经典四阶龙格库塔法解一阶微分方程组.doc

数值计算课程设计1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析（1-1）， （1-2） （1-3） （1-4）（1-5）（1-6）（1-7）（1-8）（1-9） （1-10）经过循环计算由 推得 每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。1.2经典四阶龙格库塔法解一阶微分方程流程图 图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码： #include #include using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial3, double resu3,double h) double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1; t0=initial0;x0=initial1;y0=initial2; f1=f(t0,x0,y0); g1=g(t0,x0,y0); f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2); f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2, x0+h*f2/2,y0+h*g2/2); f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6;resu0=t0+h;resu1=x1;resu2=y1;int main() double f(double t,double x, double y); double g(double t,double x, double y); double initial3,resu3; double a,b,H; double t,step; int i; coutinitial0initial1initial2; coutab; coutstep; coutsetiosflags(ios:right)setiosflags(ios:fixed)setprecision(10); H=(b-a)/step; cout initial0setw(18)initial1setw(18)initial2endl; for(i=0;istep;i+) RK4( f,g ,initial, resu,H); coutresu0setw(20)resu1setw(20)resu2endl; initial0=resu0;initial1=resu1;initial2=resu2; return(0);double f(double t,double x, double y) double dx; dx=x+2*y;return(dx);double g(double t,double x, double y)double dy;dy=3*x+2*y;return(dy);1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示： 应用所编写程序计算所给例题： 其中初值为求解区间为0,0.2。 计算结果为: 图1-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图452.高斯列主元法解线性方程组2.1高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程：for k=1 to n-1 dofor i=k+1 to nl=a(i,k)/a(k,k)for j=k to n doa(i,j)=a(i,j)-l*a(k,j) end %end of for jb(i)=b(i)-l*b(k) end %end of for iend %end of for k最后得到A，b可以构成上三角线性方程组接着使用回代法求解上三角线性方程组 因为高斯消元要求a（k,k）0（k=1,2,3n-1）这就需要对高斯消元过程进行完善，即使用高斯列主元法：其步骤为： 找主元：计算 ,并记录其所在行r , 交换第r行与第k行； 以第k行为工具行处理以下各行，使得从第k列的第k+1行到第n行的元素全部为0； 得到增广矩阵的上三角线性方程组；使用回代法对上三角线性方程组进行求解2.2高斯列主元法解线性方程组流程图FFTT开始输入A交换A中r、k两行输出x结束图2-1 高斯列主元法解线性方程组流程图2.3高斯列主元法解线性方程组程序代码#include#include #define N 3using namespace std;void main()int i,j,k,n,p;float t,s,m,aNN,bN,xN;cout请输入方程组的系数endl;for(i=0;iN;i+)for(j=0;jaij;cout请输入方程组右端的常数项:endl;for(i=0;ibi;for(j=0;jN-1;j+) p=j; for(i=j+1;ifabs(apj) p=i; if(p!=j) /交换第p行与第j行for(k=0;kN;k+)t=ajk;ajk=apk;apk=t; /交换常数项的第p行与第j行 t=bp; bp=bj; bj=t; /把系数矩阵第j列ajj下面的元素变为0for(i=j+1;iN;i+) m=-aij/ajj; for(n=j;n=0;i-) s=0.0;for(j=N-1;ji;j-)s=s+aij*xj; xi=(bi-s)/aii; cout方程组的解如下:endl;for(i=0;i=N-1;i+) coutxiendl;2.4高斯列主元法解线性方程组程序调试结果图示：求解下列方程组 图2-2 高斯列主元法解线性方程组程序算法调试图 3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明 牛顿法主要思想是用 进行迭代 。因此首先需要算出的雅可比矩阵，再求过求出它的逆，当它达到某个精度时即停止迭代。 具体算法如下：首先设已知。：计算函数 （3-1） ：计算雅可比矩阵 （3-2） （3-3）：求线性方程组的解。：计算下一点重复上述过程。3.2牛顿法解非线性方程组流程图图3-1 牛顿法解非线性方程组流程图3.3牛顿法解非线性方程组程序代码#include#include#define N 2 / 非线性方程组中方程个数、未知量个数 #define Epsilon 0.0001 / 差向量1范数的上限#define Max 100 /最大迭代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);/计算向量函数的因变量向量yyNvoid ffjacobian(float xxN,float yyNN);/计算雅克比矩阵yyNNvoid inv_jacobian(float yyNN,float invNN);/计算雅克比矩阵的逆矩阵invvoid newdundiedai(float x0N, float invNN,float y0N,float x1N);/由近似解向量 x0 计算近似解向量 x1float x0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;/如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0读入初始近似解向量/for( i=0;ix0i;cout初始近似解向量：endl;for (i=0;iN;i+)coutx0i ; coutendl;coutendl;doiter=iter+1;cout第 iter 次迭代开始endl;/计算向量函数的因变量向量 y0ff(x0,y0);/计算雅克比矩阵 jacobianffjacobian(x0,jacobian);/计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian(jacobian,invjacobian);/由近似解向量 x0 计算近似解向量 x1newdundiedai(x0, invjacobian,y0,x1);/计算差向量的1范数errornormerrornorm=0;for (i=0;iN;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornormEpsilon) break;for (i=0;iN;i+)x0i=x1i; while (iterMax);return 0;void ff(float xxN,float yyN)float x,y; int i; x=xx0;y=xx1;yy0=x*x-2*x-y+0.5;yy1=x*x+4*y*y-4; cout向量函数的因变量向量是： endl;for( i=0;iN;i+) coutyyi ; coutendl;coutendl;void ffjacobian(float xxN,float yyNN) float x,y; int i,j;x=xx0;y=xx1;/jacobian have n*n elementyy00=2*x-2;yy01=-1;yy10=2*x;yy11=8*y; cout雅克比矩阵是： endl;for( i=0;iN;i+) for(j=0;jN;j+) coutyyij ;coutendl; coutendl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k; cout开始计算雅克比矩阵的逆矩阵 ：endl; for (i=0;iN;i+) for(j=0;jN;j+) augij=yyij; for(j=N;jN2;j+)if(j=i+N) augij=1;else augij=0; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl;coutendl;for (i=0;iN;i+) for (k=i+1;kN;k+) L=-augki/augii;for(j=i;jN2;j+) augkj=augkj+L*augij; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl;cout0;i-) for (k=i-1;k=0;k-) L=-augki/augii;for(j=N2-1;j=0;j-) augkj=augkj+L*augij; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl;cout=0;i-) for(j=N2-1;j=0;j-)augij=augij/augii;for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl; for(j=N;jN2;j+) invij-N=augij;coutendl;cout雅克比矩阵的逆矩阵： endl;for (i=0;iN;i+) for(j=0;jN;j+) coutinvij ; coutendl;coutendl;void newdundiedai(float x0N, float invNN,float y0N,float x1N)int i,j;float sum=0;for(i=0;iN;i+) sum=0; for(j=0;jN;j+) sum=sum+invij*y0j;x1i=x0i-sum; cout近似解向量：endl;for (i=0;iN;i+) coutx1i ; coutendl;cout=k 的近似积分结果逼近表，并以R(j+1,j+1)为最终解来逼近积分。R(j,0)=T(j), j=0,T(j)为区间逐次减半递推梯形求积分公式算出的结果；R(j,1)=S(j), j=1,S(j)为区间逐次减半递推辛普森求积分公式算出的结果；（4-1）R(j,2)=B(j), j=2,B(j)为递推布尔求积分公式算出的结果； （4-2）生成的逼近表，并以为最终解来逼近积分（4-3）逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为当时，程序在第行结束。4.2龙贝格积分算法流程图 图4-1 龙贝格积分算法流程图4.3龙贝格积分算法程序代码#includeusing namespace std;#include#include#define f(x) (sin(x) /列举函数#define N_H 20#define MAX 10 /最大迭代次数#define a 0 /所求积分的上下限#define b 1#define epsilon 0.0001 /所需精度double Romberg(double aa,double bb,long int n) int i; double sum,h=(bb-aa)/n;sum=0; for(i=1;in;i+) sum+=f(aa+i*h); sum+=(f(aa)+f(bb)/2; return (h*sum);void main() int i; long int n=N_H,m=0; double T2MAX+1; T10=Romberg(a,b,n); n*=2; for(m=1;mMAX;m+) for(i=0;im;i+) T0i=T1i; T10=Romberg(a,b,n); n=n*2; for (i=1;i=m;i+) T1i=T1i-1+(T1i-1-T0i-1)/(pow(2,2*m)-1); if(T0m-1-T1m)epsilon) coutT=T1mendl;return; 4.4龙贝格积分算法调试图求在区间上的精度为0.0001的积分 图4-2 龙贝格积分程序算法调试图5.三次样条插值算法5.1三次样条插值基本算法说明： 表5-1三次样条插值基本算法说明策略描述包含和的方程(i)三次紧压样条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”），(iii)外挂到端点若已知N+1个点的及其一阶导数的边界条件S(a)=和S(b)=,则存在唯一的三次样条曲线。构造并求解下列线性方程组（5-1） （5-2）当得到系数 后，可利用如下公式计算样条函数 （5-3） 为了更有效的计算，每个三次多项式 可表示成嵌套乘的形式，也可以写为如下形式：（5-4）5.2三次样条插值算法程序代码#include #include #define MAX 4 double *diff(double X,int n) int i=0; double *H=NULL; H=(double*)malloc(n-1)*sizeof(double); for(i=1;i=n-1;i+) Hi-1=Xi-Xi-1; return H; double *divide(double Y,int N,double H) int i=0; double *D=NULL; D=(double*)malloc(N*sizeof(double); for(i=0;iN;i+) Di=Yi/Hi; return D; main() doubleXMAX=0,1,2,3,YMAX=0,0.5,2.0,1.5,SMAXMAX=0,temp=0.0,MMAX=0; int N=MAX-1,i=0,k=0; double AMAX-1-2=0,BMAX-1-1=0,CMAX-1-1=0; double dx0=0.2,dxn=1.0; double *H=NULL,*D=NULL,*U=NULL; printf(求解过点（0，0.0），（1，0.5），（2，2.0）和（3，1.5）n而且一阶导数边界条件S(0)=0.2和S(2)=-1的三次压紧样条曲线nn); H=diff(X,MAX); D=divide(diff(Y,MAX),N,H); for(i=1;iN-2;i+) Ai=Hi+1; for(i=0;iN-1;i+) Bi=2*(Hi+Hi+1); for(i=1;iN-1;i+) Ci=Hi+1; U=diff(D,N); for(i=0;iN;i+) Ui=Ui*6; B0=B0-H0/2; U0=U0-3*(D0-dx0); BN-2=BN-2-HN-1/2; UN-2=UN-2-3*(dxn-DN-1); for(k=2;k=1;k-) Mk=(Uk-1-Ck-1*Mk+1)/Bk-1; M0=3*(D0-dx0)/H0-M0/2; MN=3*(dxn-DN-1)/HN-1-MN-1/2; for(k=0;k=N-1;k+) Sk0=(Mk+1-Mk)/(6*Hk); Sk1=Mk/2; Sk2=Dk-Hk*(2*Mk+Mk+1)/6; Sk3=Yk; printf(求得的三次紧压样条曲线的矩阵S为:n); for(i=0;iMAX-1;i+) for(k=0;kMAX;k+) printf(%lft,Sik); printf(n); for(i=0;i x1=0:.01:1; y1=polyval(S(1,:),x1-X(1);x2=1:.01:2; y2=polyval(S(2,:),x2-X(2);x3=2:.01:3; y3=polyval(S(2,:),x3-X(3);plot(x1,y1,x2,y2,x3,y3,X,Y,.)运行结果如下： 图5-2 三次样条插值作图 6.M次多项式曲线拟合6.1 M次多项式曲线拟合算法说明：设有N个点，横坐标是确定的。最小二乘法抛物线系数表示为，（6-1）求解A，B和C的线性方程组为（6-2）整体上考虑近似函数同所给数据点(i=0,1,m)误差 (i=0,1,m)的大小，常用的方法有以下三种：一是误差(i=0,1,m)绝对值的最大值，即误差向量的范数；二是误差绝对值的和，即误差向量r的1范数；三是误差平方和的算术平方根，即误差向量r的2范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2范数的平方，因此在曲线拟合中常采用误差平方和来 度量误差 (i=0，1，m)的整体大小。数据拟合的具体作法是：对给定数据 (i=0,1,，m)，在取定的函数类中,求,使误差(i=0,1,m)的平方和最小。6.2 M次多项式曲线拟合流程图FT开始输入x0、输出极小值点和极小值结束图6-1 M次曲线多项式曲线拟合流程图6.3 M次多项式曲线拟合算法程序代码#include #include #define MAX 20 void inv(double XMAXMAX,int n,double EMAXMAX) int i=0,j=0,k=0; double temp=0.0; for(i=0;iMAX;i+) for(j=0;jMAX;j+) if(i=j) Eij=1; for(i=0;in-1;i+) temp=Xii; for(j=0;jn;j+) Xij=Xij/temp; Eij=Eij/temp; for(k=0;kn;k+) if(k=i) continue; temp=-Xii*Xki; for(j=0;jn;j+) Xkj=Xkj+temp*Xij; Ekj=Ekj+temp*Eij; main() int n=0,M=0,i=0,j=0,k=0; double XMAX=0,YMAX=0,FMAXMAX=0,BMAX=0; double AMAXMAX=0,BFMAXMAX=0,EMAXMAX=0,CMAX=0; printf(tttM次多项式曲线拟合nn输入拟合的点的个数:); scanf(%d,&n); printf(n输入%d个点的X坐标序列:n,n); for(i=0;in;i+) scanf(%lf,&Xi); printf(n输入%d个点的Y坐标序列:n,n); for(i=0;in;i+) scanf(%lf,&Yi); printf(n输入拟合次数:); scanf(%d,&M); for(i=0;in;i+) for(k=1;k=M+1;k+) Fik-1=pow(Xi,k-1); for(i=0;in;i+) for(j=0;jM+1;j+) BFji=Fij; for(i=0;iM+1;i+) for(j=0;jM+1;j+) for(k=0;kn;k+) Aij+=BFik*Fkj; for(i=0;iM+1;i+) for(j=0;jn;j+) Bi+=BFij*Yj; inv(A,n,E); for(i=0;iM+1;i+) for(j=0;j=0;i-) printf(%lft,Ci); printf(n拟合后的%d次多项式为:n,M); printf(nP(x)=); for(i=M;i=0;i-) if(i=0) printf(%+.3lfn,Ci); else printf(%+.3lf*x%d,Ci,i); getchar(); 6.4 M次多项式曲线拟合算法程序调试结果：求解（-3，3），（0，1），（2，1），（4，3）四个点拟合2次后的2后的多项式图6-2 M次曲线多项式曲线拟合程序算法调试图6.5 M次多项式曲线拟合x=-3:0:2:4;y=3:1:1:3;plot(x,y,b)图6-3 M次多项式曲线拟合画图7.二分法解非线性方程7.1 二分法解非线性方程算法说明：.是起始区间，是中点。.是第二个区间，它包含零点r，同时是中点，区间的宽度范围是的一半。.得到第n个区间（包含r，并有中点）后，可构造出，它也包括r，宽度范围是的一半。 7.2二分法解非线性方程程序代码#includeusing namespace std;#define eps 0.00000000000001double f(double x) return 2-x*x;int main()double ga,gb,gc,a,b,c; cout *endl;cout 用二分法寻找方程2-x*x=0的根endl;cout 求根区间为(1,2)endl; cout *eps)c=(a+b)/2;gc=f(c);if(gc=0) break;else if(gc*gb0)a=c;ga=gc;elseb=c;gb=gc;cout 方程的根为:X=cendl; cout *endl;return 0;7.3二分法解非线性方程算法调试结果图示： 用二分法求 在1,2区间上的根图7-2二分法解非线性方程程序算法调试图8.不动点法解非线性方程8.1 不动点法解非线性方程算法概述将改写成对进行迭代，即 其中判断是否成立，成立则返回，不成立就重复以上步骤8.2 不动点法说明函数g(x)的不动点是指一个实数P，满足P=g(P)。从图形角度分析，函数y=g(x)的不动点是y=g(x)和y=x的交点。求不动点的过程可描述如下：第一步：将x0带入迭代公式，作为第一次迭代结果x1第二步：随着不断的迭代，迭代数值会越来越接近不动点的值x0，程序中变量 的类型小数点后的位数是一定的，所以，随着不断的迭代，会出现在误 差线内的两数，那么，此时的可以近似看做方程的根。8.3不动点法解非线性方程程序代码：#includeusing namespace std;#include#define MAX 20#define eps 2.2204e-16double g(double x)return 1-x*x/2+x;main()double PMAX=0,err=0.0,relerr=0.0,tol=0.0,p=0.0,p0=0.0;int k=0,max1=0,i=0;cout不动点法解非线性方程f(x)=1-x2/2+xendl;cout方程在0，1上有解，初始值为p0=0endl;P0=p0=0;max1=100;tol=0.001; for(k=2;k=max1;k+) Pk-1=g(Pk-2);err=fabs(Pk-1-Pk-2);relerr=err/(fabs(Pk-1+eps);p=Pk-1;if(errtol)|(relerrtol)break;if(k=max1) cout迭代次数超过允许的最大迭代次数！endl;cout不动点的近似值为:pendl;cout程序迭代次数为kendl;cout近似值的误差为:errendl;cout求解不动点近似值的序列:endl;for(i=0;ik;i+) cout Pi;coutendl;8.4不动点法解非线性方程程序调试图示用不动点法求解图8-1不动点法解非线性方程程序算法调试图9.牛顿-拉弗森迭代解非线性方程 9.1牛顿-拉弗森迭代解非线性方程算法说明： 存在数pa,b,满足f(p)=0.如果存在一个数对任意近似值，使得由如下迭代定义的序列 收敛到p： （9-1）其中并被称为牛顿-拉弗森迭代函数。由于f(p)=0，显然g(p)=p.这样通过寻找函数g(x)的不动点，可以实现寻找方程f(x)=0的根牛顿-拉弗森迭代。9.2 牛顿-拉弗森迭代解非线性方程流程图NY图9-1 牛顿-拉弗森迭代解非线性方程流程图9.3牛顿-拉弗森迭代解非线性方程程序代码：#include#includeusing namespace std;double f(double x)return x*x*x-3*x*x-x+1; double df(double x)return 3*x*x-6*x-1;int main() int k=0,max1=0; double p0,delta,epsilon,p1,err,relerr,y;cout牛顿拉弗森法解非线性方程f(x)=x2-1endl;cout初始值为p0=4.8endl; p0=4.8,delta=0.0001,epsilon=0.0001,max1=100;for(k=1;k=max1;k+)p1=p0-f(p0)/df(p0);err=fabs(p1-p0);relerr=2*err/(fabs(p1)+delta);p0=p1;y=f(p0);if(errdelta)|(relerrdelta)|(fabs(y)epsilon)break;cout方程的根的近似值为:p0endl;cout方程的根的误差估计为:errendl;cout迭代次数为:kendl;cout方程在p0点的函数值f(p0):yendl;return 0;9.4牛顿-拉弗森迭代解非线性方程算法程序调试：用牛顿-拉弗森迭代解非线性方程图9-2牛顿-拉弗森迭代解非线性方程程序算法调试图10.拉格朗日插值10.1 拉格朗日插值算法说明有n+1个点的次数最高为n项的多项式的构造方法,其中拉格朗日（10-1）的系数多项式对每个固定的k，拉格朗日系数多项式具有性质（10-2）10.2 拉格朗日插值算法程序#include using namespace std; double func(double X,int k,double x,int n); int main() double Sn=0; int n; coutn; double*x=(double*)malloc(n*sizeof(double); double*y=(double*)malloc(n*sizeof(double); double X; int i; for(i=0;in;i+) cout请输入xi+1,yi+1:xiyi; coutX; for(i=0;in;i+) Sn=Sn+func(X,i,x,n)*yi; cout通过拉格朗日插值公式所得的结果：当x=X时，y=Snendl; return 0; double func(double X,int k,double x,int n) int i; double Pn=1; double p; for(i=0;in;i+) if(i=k) continue; else p=(X-xi)/(xk-xi); Pn=Pn*p; return Pn; 10.3 拉格朗日插值算法程序调试：求解（1，2），（8，9），（0，0），（101）四个点的拉格朗日插值图10-1程序算法调试图设计体会及今后的改进意见这次课程设计我大约花了两个周的时间，开始的时候还觉得挺简单的，后来当我真正的开始做的时候，才发现这不是一件简单的事。然而这些算法都来自课本，对于课本理论知识本来就薄弱的我来说，这更加是难上加难不过我相信，只要认真做，就没有完成不了的。看了很久的课本再加上在网上查了很多相关的资料后，终于有了一些眉目。程序的调试阶段也是非常困难的，有时觉得程序明明就是对的，但它就是运行不出来，错误提示也很不准确，在这些问题当中，主要以程序的逻辑问题，语法问题等问题为主。只有把程序读懂，一句一句的看，观察每个语句的运行结果，才能保证程序的每一步计算都是正确的。通过这次典型数值算法的C+语言程序设计，我知道了很多以前在C+课本上没有学到的知识，知道了遇到问题不能自暴自弃，要冷静客观的分析问题出现在哪里，从而更好的解决问题。更重要的是，要学会与