数值计算课程设计经典四阶龙格库塔法解一阶微分方程.doc

《数值计算课程设计经典四阶龙格库塔法解一阶微分方程.doc》由会员分享，可在线阅读，更多相关《数值计算课程设计经典四阶龙格库塔法解一阶微分方程.doc（40页珍藏版）》请在沃文网上搜索。

1、数值计算课程设计1、 经典四阶龙格库塔法解一阶微分方程1.1、 算法说明龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。这种算法可以描述为，自初始点开始，利用下面的计算方法生成近似序列(1-1)1.2、 经典四阶龙格库塔法解一阶微分方程算法流程图图1-1 经典四阶龙格库塔法解一阶微分方程算法流程图1.3、 经典四阶龙格库塔法解一阶微分方程程序调试图1-2 经典四阶龙格库塔法解一阶微分方程程序调试1.4、

2、经典四阶龙格库塔法解一阶微分方程代码#include #include using namespace std;/f为函数的入口地址，x0、y0为初值，xn为所求点，step为计算次数double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, int step ) double k1,k2,k3,k4,result; double h=(xn-x0)/step; if(step=0) return(y0); if(step=1) k1=f(x0,y0); k2=f(x0+h/2, y0

3、+h*k1/2); k3=f(x0+h/2, y0+h*k2/2); k4=f(x0+h, y0+h*k3); result=y0+h*(k1+2*k2+2*k3+k4)/6; else double x1,y1; x1=xn-h; y1=Runge_Kuta(f, x0, y0, xn-h,step-1); k1=f(x1,y1); k2=f(x1+h/2, y1+h*k1/2); k3=f(x1+h/2, y1+h*k2/2); k4=f(x1+h, y1+h*k3); result=y1+h*(k1+2*k2+2*k3+k4)/6; return(result);int main()do

4、uble f(double x, double y); double x0,y0; double a,b;/ int step; coutx0y0; coutab; /double x0=0,y0=1; double x,y,step; int i; coutstep; /step=0.1; cout.precision(10); for(i=0;i=(b-a)/step;i+)x=x0+i*step; coutsetw(8)xsetw(18)Runge_Kuta(f,x0,y0,x,i)endl; return(0); double f(double x, double y) double

5、r; r=(x-y)/2; return(r); 2、 高斯列主元法解线性方程组2.1、 算法说明首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。对第元素，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元法将新得到的消去第i行以下的元素。一次进行直到。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。2.2、 高斯列主元算法流程图i=0;j=0NY开始输入未知数个数mim*(m+1)max j=i:max=|*(head+(m+1)*i+i)|;kmaxmax =|*(head+(m+1)*k+i)|;max i=

6、k;NYYKmK=k+1;YNmax i!=itemp=*(head*(m+1)*i+k);*(head*(m+1)*i+k)=*(head*(m+1)*max i+k)*(head*(m+1)*max i+k)=tempk=0N- 39 -k=m+1Y NYi!=jk=0*(head*(m+1)*i+k)=*(head+(m+1)*j+k)*(head*(m+1)*j+i) *(head*(m+1)* i+k)/(*head+(m+1)*i+i)k=m+1YNj=j+1YjmYNi=0NY结束*(head+(m+1)*i+m)=*(head+(m+1)*i+m)/(*head+(m+1)*i

7、+m)；i=i+1；图2-1 算法流程图输出*(head+(m+1)*i+m)jm2.3、 高斯列主元程序调试对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3，输入3行5列的增广矩阵,运行界面为：图2-2 高斯列主元程序调试2.4、 高斯列主元算法代码#include #includeusing namespace std;void load();const N=20;float aNN;int m;int main()int i,j;int c,k,n,p,r;float xN,lNN,s,d;coutm;coutendl;cout请按顺序输

8、入增广矩阵a:endl;load();for(i=0;im;i+) for(j=i;jfabs(aii)?j:i; /*找列最大元素*/for(n=0;nm+1;n+) s=ain; ain=acn; acn=s; /*将列最大数防在对角线上*/for(p=0;pm+1;p+)coutaipt;coutendl;for(k=i+1;km;k+) lki=aki/aii; for(r=i;r=0;i-) d=0;for(j=i+1;jm;j+)d=d+aij*xj;xi=(aim-d)/aii;cout该方程组的解为:endl;for(i=0;im;i+)coutxi=xit; /system(

9、pause); return 0; void load()int i,j;for(i=0;im;i+)for(j=0;jaij;3、 牛顿法解非线性方程组3.1、 算法说明设已知。第1步：计算函数 （3-1）第2步：计算雅可比矩阵 （3-2）第3步：求线性方程组 的解。第4步：计算下一点 （3-3）重复上述过程。3.2、 牛顿法解非线性方程组算法流程图图3-1 算法流程图3.3、 牛顿法解非线性方程组算法程序调试图3-2牛顿法解非线性方程组算法程序调试应用本程序解方程组， 初始近似值x0,y0分别为2.00和0.25，经过3次迭代求出X(1)=1.900691和X(2)=0.311213。图3

10、-2牛顿法解非线性方程组算法程序运行结果3.4、 牛顿法解非线性方程组算法程序代码#include#include#define N 2 / 非线性方程组中方程个数、未知量个数 #define epsilon 0.0001 / 差向量1范数的上限#define max 10 /最大迭代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);void ffjacobian(float xxN,float yyNN);void inv_jacobian(float yyNN,float invNN)

11、;void newdim(float x0N, float invNN,float y0N,float x1N);float x0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,iter=0;cout初始解向量：endl;for (i=0;iN;i+) coutx0i ; coutendl;do iter=iter+1;cout第 iter 次迭代开始：endl;/jis uan xiang liang han shu zhi-yin bian liang xiang liang y0ff(x0,y0);/jis uan

12、jacobian ju zhen jacobianffjacobian(x0,jacobian);/jis uan jacobian ju zhen de ni juzhen invjacobianinv_jacobian(jacobian,invjacobian);/you jie xiang liang x0 ji suan jie xiang liang x1newdim(x0, invjacobian,y0,x1);/ji suan cha xiang liang de 1 fan shuerrornorm=0;for (i=0;iN;i+) errornorm=errornorm+f

13、abs(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;in

14、t i,j;x=xx0;y=xx1;yy00=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

15、 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

16、;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 newdim(float x0N, float invNN,flo

17、at 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 endl; 4、 龙贝格求积分算法4.1、 算法说明生成的逼近表，并以为最终解来逼近积分 （4-1）逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为 （4-2）当时，程序在第行结束。4.2、 龙贝格求积分

18、算法流程图图4-1 算法流程图图4-1 算法流程图4.3、 龙贝格求积分算法程序调试我们以求解积分，精度为0.0001，最高迭代10次为例，对所编写的龙贝格求积分算法程序进行编译和链接，经执行后得如下所示的窗口图4-2龙贝格求积分算法程序调试说明：应用Romberg算法求在区间上的精度为0.0001的积分为0.494508。4.4、 龙贝格求积分算法代码#include#includeusing namespace std;#define f(x) sin(x*x) /举例函数#define epsilon 0.0001 /精度#define MAXREPT 10 /迭代次数,到最后仍达不到精

19、度要求,则输出T(m=10).double Romberg(double aa, double bb) /aa,bb 积分上下限 int m, n;/m控制迭代次数, 而n控制复化梯形积分的分点数. n=2m double h, x; double s, q; double ep; /精度要求 double *y = new doubleMAXREPT;/为节省空间,只需一维数组 /每次循环依次存储Romberg计算表的每行元素,以供计算下一行,算完后更新double p;/p总是指示待计算元素的前一个元素(同一行) /迭代初值 h = bb - aa; y0 = h*(f(aa) + f(b

20、b)/2.0; m = 1; n = 1; ep = epsilon + 1.0; /迭代计算 while (ep = epsilon) & (m MAXREPT) /复化积分公式求T2n(Romberg计算表中的第一列),n初始为1,以后倍增 p = 0.0; for (int i=0; in; i+)/求Hn x = aa + (i+0.5)*h; p = p + f(x); p = (y0 + h*p)/2.0;/求T2n = 1/2(Tn+Hn),用p指示 /求第m行元素,根据Romberg计算表本行的前一个元素(p指示), /和上一行左上角元素(yk-1指示)求得. s = 1.0;

21、 for (int k=1; k=m; k+) s = 4.0*s; q = (s*p - yk-1)/(s - 1.0); yk-1 = p; p = q; p = fabs(q - ym-1); m = m + 1; ym-1 = q; n = n + n; h = h/2.0; return (q);int main() double a,b; coutRomberg积分,请输入积分范围a,b:ab; cout积分结果:Romberg(a, b)endl; return 0;5、 三次样条插值算法5.1 三次样条插值算法说明表表5-1 三次样条插值算法说明表策略描述包含和的方程(i)三次

22、紧压样条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”），(iii)外挂到端点(iv) 是靠近端点的常量，(v)在每个端点处指定，5.2、 三次样条插值算法（压紧样条）程序调试我们将所编写的程序三次样条插值算法（压紧样条）程序进行调试图5-1三次样条插值算法（压紧样条）程序输入界面、运行结果图5-2三次样条插值算法程序运行结果（a） 图5-2三次样条插值算法程序运行结果 （b）运行结果分析： （5-1）作图程序（Matlab）：x1=0:0.01:1;y1=0.06*(x1 - 1).3 + 0.42*(x1 - 0).3 - 0.06*(x1 -

23、1) + 0.08*(x1 - 0);x2=1:0.01:2;y2=-0.42*(x2 - 2).3 - 0.62*(x2 - 1).3 - 0.08*(x2 - 2) + 2.62*(x2 - 1);x3=2:0.01:3;y3=0.62*(x3 - 3).3 + 0.06*(x3 - 2).3 - 2.62*(x3 - 3) + 1.44*(x3 - 2);X=0 1 2 3;Y=0 0.5 2 1.5;plot(x1,y1,x2,y2,x3,y3,X,Y,*)gtext(S1)gtext(S2)gtext(S3)图形为：图5-3 三次样条插值算法Matlab作图分析5.3、 三次样条插值

24、算法（压紧样条）代码#include#includeusing namespace std;const int max = 50;float xmax, ymax, hmax;float cmax, amax, fxymmax;float f(int x1, int x2, int x3) float a = (yx3 - yx2) / (xx3 - xx2); float b = (yx2 - yx1) / (xx2 - xx1); return (a - b)/(xx3 - xx1); /求差分void cal_m(int n) /用追赶法求解出弯矩向量M float Bmax; B0 =

25、 c0 / 2; for(int i = 1; i n; i+) Bi = ci / (2 - ai*Bi-1); fxym0 = fxym0 / 2; for(i = 1; i = 0; i-) fxymi = fxymi - Bi*fxymi+1;void printout(int n);int main() int n,i; char ch; do coutn; for(i = 0; i = n; i+) coutPlease put in Xixi; coutPlease put in Yiyi; for(i = 0; i n; i+) /求步长 hi = xi+1 - xi; cou

26、tt; switch(t) case 1:cout输入 Y0 Ynf0f1; c0 = 1; an = 1; fxym0 = 6*(y1 - y0) / (x1 - x0) - f0) / h0; fxymn = 6*(f1 - (yn - yn-1) / (xn - xn-1) / hn-1; break; case 2:cout输入 Y0 Ynf0f1; c0 = an = 0; fxym0 = 2*f0; fxymn = 2*f1; break; default:cout不可用n;/待定 ; for(i = 1; i n; i+) fxymi = 6 * f(i-1, i, i+1);

27、for(i = 1; i n; i+) ai = hi-1 / (hi + hi-1); ci = 1 - ai; an = hn-1 / (hn-1 + hn); cal_m(n); coutn输出三次样条插值函数：n; printout(n); coutch; while(ch = y | ch = Y); return 0;void printout(int n)coutsetprecision(6); for(int i = 0; i n; i+)couti+1: xi , xi+1nt;coutSi+1 0) cout-t*(x - xi+1)3; else cout-t*(x -

28、xi+1 0) cout + t*(x - xi)3; else cout - t*(x - xi)3; cout 0) cout- t*(x - xi+1); else cout- -t*(x - xi+1 0) cout + t*(x - xi); else cout - -t*(x - xi); coutendl; coutendl;6、M次多项式曲线拟合6.1、算法说明设有N个点，横坐标是确定的。最小二乘抛物线的系数表示为 （6-1）求解A，B和C的线性方程组为 （6-2）6.2、 M次多项式曲线拟合算法流程图图6-1 算法流程图6.3、 M次多项式曲线拟合算法程序调试我们按命令依次输

29、入命令如下命令后，得程序执行结果如下图6-2 算法调试图作图程序： 图形为：X=-3 -1 1 3;Y=15 5 1 5;x=-3:0.01:3;y=0.875*x.2-1.7*x+2.145;plot(X,Y,go,x,y,r-)gtext(拟合曲线)图6-3 图形分析6.4、 M次多项式曲线拟合算法代码#include#include#define MAX 20using namespace std;/求解任意可逆矩阵的逆，X为待求解矩阵，E为全零矩阵，非单位矩阵，也可以是单位矩阵void inv(double XMAXMAX,int n,double EMAXMAX)int i,j,k;

30、double temp=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;int main()int n,M,i,j,k;double XMAX=0,YMAX=0,FMAXMAX=0,BMAX=0;double AMAXMAX=0,BFMAX

31、MAX=0,CMAX=0,EMAXMAX=0;coutn;coutn请输入n个点的X坐标序列:n;for(i=0;iXi;coutn请输入n个点的Y坐标序列:n;for(i=0;iYi;coutM;for(i=0;in;i+)for(k=1;k=M+1;k+)Fik-1=pow(Xi,k-1);/求F的转置for(i=0;in;i+)for(j=0;jM+1;j+)BFji=Fij; /计算其转置的BF与F的乘 for(i=0;iM+1;i+)for(j=0;jM+1;j+)for(k=0;kn;k+)Aij+=BFik*Fkj;/计算F的转置BF与Y的乘for(i=0;iM+1;i+)for(j=0;jn;j+)Bi+=BFij*Yj;/调用inv函数求解矩阵A的逆矩阵Einv(A,n,E);/计算A的逆BF与B的乘for(i=0;iM+1;i+)for(j=0;jn;j+)Ci+=Eij*Bj;coutn拟合后的M=0;i-)coutCit;coutn拟合后的M次多项式为:n;cout=0;i-)