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

数值计算课程设计11.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析),,(1kkkyxtff，)2,2,2(112ghyfhxhtffkkk)2,2,2(223ghyfhxhtffkkk),,(334hgyhfxhtffkkk),,(1kkkyxtgg)2,2,2(112ghyfhxhtggkkk)2,2,2(223ghyfhxhtggkkk),,(334hgyhfxhtggkkk)22(6)22(64321143211gggghyyffffhxxkkkk1kktth经过循环计算由推得……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为NOh,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。000,,txy111222,,,,txytxy（1-1）（1-2）（1-3）（1-4）（1-5）（1-6）（1-7）（1-8）（1-9）（1-10）经典四阶龙格库塔法解一阶微分方程组21.2经典四阶龙格库塔法解一阶微分方程流程图图1-1经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码：#includeiostream#includeiomanipusingnamespacestd;voidRK4(double(*f)(doublet,doublex,doubley),double(*g)(doublet,doublex,doubley),doubleinitial[3],doubleresu[3],doubleh){doublef1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial[0];x0=initial[1];y0=initial[2];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,数值计算课程设计3x0+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;resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}intmain(){doublef(doublet,doublex,doubley);doubleg(doublet,doublex,doubley);doubleinitial[3],resu[3];doublea,b,H;doublet,step;inti;cout输入所求微分方程组的初值t0,x0,y0:;cininitial[0]initial[1]initial[2];cout输入所求微分方程组的微分区间[a,b]:;cinab;cout输入所求微分方程组所分解子区间的个数step:;cinstep;coutsetiosflags(ios::right)setiosflags(ios::fixed)setprecision(10);H=(b-a)/step;coutinitial[0]setw(18)initial[1]setw(18)initial[2]endl;for(i=0;istep;i++){RK4(f,g,initial,resu,H);coutresu[0]setw(20)resu[1]setw(20)resu[2]endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}doublef(doublet,doublex,doubley){doubledx;经典四阶龙格库塔法解一阶微分方程组4dx=x+2*y;return(dx);}doubleg(doublet,doublex,doubley){doubledy;dy=3*x+2*y;return(dy);}1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示：应用所编写程序计算所给例题：其中初值为求解区间为[0,0.2]。计算结果为:图1-2经典四阶龙格库塔法解一阶微分方程算法程序调试图yxdtdyyxdtdx2324)0(6)0(yx数值计算课程设计52.高斯列主元法解线性方程组2.1高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程：fork=1ton-1dofori=k+1tonl=a(i,k)/a(k,k)forj=ktondoa(i,j)=a(i,j)-l*a(k,j)end%endofforjb(i)=b(i)-l*b(k)end%endofforiend%endoffork最后得到A，b可以构成上三角线性方程组接着使用回代法求解上三角线性方程组因为高斯消元要求a（k,k）≠0（k=1,2,3……n-1）这就需要对高斯消元过程进行完善，即使用高斯列主元法：其步骤为：①找主元：计算1max||()kpkapkn,并记录其所在行r,||rka1max||()kpkapkn②交换第r行与第k行；③以第k行为工具行处理以下各行，使得从第k列的第k+1行到第n行的元素全部为0；④得到增广矩阵的上三角线性方程组；⑤使用回代法对上三角线性方程组进行求解2.2高斯列主元法解线性方程组流程图高斯列主元法解线性方程组6图2-1高斯列主元法解线性方程组流程图FFTT开始输入Ak1rkikkiaamax?0rka交换A中r、k两行ijkkkjikijaaaaankki,,2,11,,2,1nkkj?1nkkkknkjjkjknxaxaa111,,1,nnk输出x结束kk1数值计算课程设计72.3高斯列主元法解线性方程组程序代码#includeiostream#includecmath#defineN3usingnamespacestd;voidmain(){inti,j,k,n,p;floatt,s,m,a[N][N],b[N],x[N];cout请输入方程组的系数endl;for(i=0;iN;i++){for(j=0;jN;j++)cina[i][j];}cout请输入方程组右端的常数项:endl;for(i=0;iN;i++)cinb[i];for(j=0;jN-1;j++){p=j;for(i=j+1;iN;i++)//寻找系数矩阵每列的最大值{if(fabs(a[i][j])fabs(a[p][j]))p=i;}if(p!=j)//交换第p行与第j行{for(k=0;kN;k++){t=a[j][k];a[j][k]=a[p][k];a[p][k]=t;}//交换常数项的第p行与第j行t=b[p];b[p]=b[j];b[j]=t;}//把系数矩阵第j列a[j][j]下面的元素变为0for(i=j+1;iN;i++){m=-a[i][j]/a[j][j];高斯列主元法解线性方程组8for(n=j;nN;n++)a[i][n]=a[i][n]+a[j][n]*m;b[i]=b[i]+b[j]*m;}}//求方程组的一个解x[N-1]=b[N-1]/a[N-1][N-1];//回代法求方程组其他解for(i=N-2;i=0;i--){s=0.0;for(j=N-1;ji;j--){s=s+a[i][j]*x[j];x[i]=(b[i]-s)/a[i][i];}}cout方程组的解如下:endl;for(i=0;i=N-1;i++)coutx[i]endl;}2.4高斯列主元法解线性方程组程序调试结果图示：求解下列方程组341-11182522321321321xxxxxxxxx数值计算课程设计9图2-2高斯列主元法解线性方程组程序算法调试图高斯列主元法解线性方程组103.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明牛顿法主要思想是用(1)()()1()()()kkkkxxFxFx(0,1,...).k进行迭代。因此首先需要算出()Fx的雅可比矩阵()Fx，再求过()Fx求出它的逆1()Fx，当它达到某个精度时即停止迭代。具体算法如下：首先设kP已知。①：计算函数12(,)()(,)kkkkkfpqFPfpq②：计算雅可比矩阵1122(,)(,)()(,)(,)kkkkkkkkkfpqfpqxyJPfpqfpqxy③：求线性方程组()()kkJPPFP的解P。④：计算下一点1kkPPP重复上述过程。（3-1）（3-2）（3-3）数值计算课程设计113.2牛顿法解非线性方程组流程图图3-1牛顿法解非线性方程组流程图3.3牛顿法解非线性方程组程序代码#includeiostream#includecmath#defineN2//非线性方程组中方程个数、未知量个数#defineEpsilon0.0001//差向量1范数的上限#defineMax100//最大迭代次数usingnamespacestd;constintN2=2*N;intmain(){voidff(floatxx[N],floatyy[N]);//计算向量函数的因变量向量yy[N]voidffjacobian(floatxx[N],floatyy[N][N]);//计算雅克比矩阵yy[N][N]voidinv_jacobian(floatyy[N][N],floatinv[N][N]);//计算雅克比矩阵的逆矩阵invvoidnewdundiedai(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]);//由近似解向量x0计算近似解向量x1floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errorno高斯列主元法解线性方程组12rm;inti,j,iter=0;//如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0读入初始近似解向量//for(i=0;iN;i++)//cinx0[i];cout初始近似解向量：endl;for(i=0;iN;i++)coutx0[i];coutendl;coutendl;do{iter=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(x1[i]-x0[i]);if(errornormEpsilon)break;for(i=0;iN;i++)x0[i