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

《经典四阶龙格库塔法解一阶微分方程组_2.docx》由会员分享，可在线阅读，更多相关《经典四阶龙格库塔法解一阶微分方程组_2.docx（10页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、经典四阶龙格库塔法解一阶微分方程组文档视界经典四阶龙格库塔法解一阶微分方程组经典四阶龙格库塔法解一阶微分方程组算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，由于它非常精准，稳定，且易于编程。1.2经典四阶龙格库塔法解一阶微分方程流程图图1-1经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码：#include#includeusingnamespacestd;voidRK4(double(*f)(doublet,doublex,doubley),double(*g)(doublet,doublex,doubley),doubleinitial

2、3,doubleresu3,doubleh)doublef1,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+

3、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;intmain()doublef(doublet,doublex,doubley);doubleg(doublet,doublex,doubley);doubleinitial3,resu3;doublea,b,H;doublet,step;inti;coutinitial0initial1initial2;coutab;coutstep;cout文档视界经典四阶龙格库塔法解一阶微分方程组

4、经典四阶龙格库塔法解一阶微分方程组end%endofforiend%endoffork最后得到A，b能够构成上三角线性方程组接着使用回代法求解上三角线性方程组由于高斯消元要求ak,k0k=1,2,3n-1这就需要对高斯消元经过进行完善，即便用高斯列主元法：其步骤为：找主元：计算,并记录其所在行r,交换第r行与第k行；以第k行为工具行处理下面各行，使得从第k列的第k+1行到第n行的元素全部为0；得到增广矩阵的上三角线性方程组；使用回代法对上三角线性方程组进行求解2.2高斯列主元法解线性方程组流程图图2-1高斯列主元法解线性方程组流程图2.3高斯列主元法解线性方程组程序代码#include#inc

5、lude#defineN3usingnamespacestd;voidmain()inti,j,k,n,p;floatt,s,m,aNN,bN,xN;coutaij;coutbi;for(j=0;jp=j;for(i=j+1;ifabs(apj)p=i;if(p!=j)/交换第p行与第j行for(k=0;k=0;i-)s=0.0;for(j=N-1;ji;j-)s=s+aij*xj;xi=(bi-s)/aii;cout2.4高斯列主元法解线性方程组程序调试结果图示：求解下列方程组图2-2高斯列主元法解线性方程组程序算法调试图3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法讲明牛顿法主要思

6、想是用进行迭代。因而首先需要算出的雅可比矩阵，再求过求出它的逆，当它到达某个精度时即停止迭代。详细算法如下：首先设已知。：计算函数3-1：计算雅可比矩阵3-23-3：求线性方程组的解。：计算下一点重复上述经过。3.2牛顿法解非线性方程组流程图图3-1牛顿法解非线性方程组流程图3.3牛顿法解非线性方程组程序代码#include#include#defineN2/非线性方程组中方程个数、未知量个数#defineEpsilon0.0001/差向量1范数的上限#defineMax100/最大迭代次数usingnamespacestd;constintN2=2*N;intmain()voidff(flo

7、atxxN,floatyyN);/计算向量函数的因变量向量yyNvoidffjacobian(floatxxN,floatyyNN);/计算雅克比矩阵yyNNvoidinv_jacobian(floatyyNN,floatinvNN);/计算雅克比矩阵的逆矩阵invvoidnewdundiedai(floatx0N,floatinvNN,floaty0N,floatx1N);/由近似解向量x0计算近似解向量x1floatx0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;inti,j,iter=0;/假如取消对x0的初始化，撤销下面两行的注释符，就能够由键盘向x0读入初始近似解向量/for(i=0;ix0i;coutvoidff(floatxxN,floatyyN)floatx,y;inti;x=xx0;y=xx1;yy0=x*x-2*x-y+0.5;yy1=x*x+4*y*y-4;coutfor(i=0;i0;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; i