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

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

1、经典四阶龙格库塔法解一阶微分方程组数值计算课程设计1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析xk1xkh(f12f22f3f4)6hyk1ykh6(g12g22g3g4)tk1tkh经过循环计算由t0,x0,y0推得t1,x1,y1t2,x2,y2?每个龙格-库塔方法都是由一个适宜的泰勒方法推导而来，使得其最终全局误差为OhN,一种折中方法是每次进行若干次函数求值，进而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，由于它非常精准，稳定，且易于编程。f1f(tk,xk,yk)，hf2f(tkh2,xkhh2f1,yk2g

2、1)hf3f(tkh2,xkhf2,yk2g2)f4f(tkh,xkhf3,ykhg3)g1g(tk,xk,yk)hg2g(tk2,xkhh2f1,yk2g1)hg3g(tk2,xkhf2h2,yk2g2)g4g(tkh,xkhf3,ykhg3)1-1)1-2)1-3)1-4)1-5)1-6)1-7)1-8)1-9)1-10)经典四阶龙格库塔法解一阶微分方程组1.2经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码：#include#includeusingnamespacestd;voidRK4(double(*f)(doublet,doublex,doub

3、ley),double(*g)(doublet,doublex,doubley),doubleinitial3,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

4、+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);图1-1经典四阶龙格库塔法解一阶微分方程流程图数值计算课程设计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;dou

5、blet,step;inti;coutinitial0initial1initial2;coutab;coutstep;cout经典四阶龙格库塔法解一阶微分方程组dx=x+2*y;return(dx);doubleg(doublet,doublex,doubley)doubledy;dy=3*x+2*y;return(dy);1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示：图1-2经典四阶龙格库塔法解一阶微分方程算法程序调试图2.高斯列主元法解线性方程组2.1高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元经过：fork=1ton-1dofori=k+1tonl数值计算课程设计由

6、于高斯消元要求a(k,k)0(k=1,2,3?n-1)这就需要对高斯消元经过进行完善，即便用高斯列主元法：其步骤为：找主元：计算max|apkk1|(pkn),并记录其所在行r,|ark|max|apkk1|(pkn)交换第r行与第k行；以第k行为工具行处理下面各行，使得从第k列的第k+1行到第n行的元素全部为0；得到增广矩阵的上三角线性方程组；使用回代法对上三角线性方程组进行求解2.2高斯列主元法解线性方程组流程图文档视界经典四阶龙格库塔法解一阶微分方程组经典四阶龙格库塔法解一阶微分方程组数值计算课程设计2.3高斯列主元法解线性方程组程序代码#include#include#defineN3

7、usingnamespacestd;voidmain()inti,j,k,n,p;floatt,s,m,aNN,bN,xN;coutaij;coutbi;for(j=0;jfabs(apj)p=i;if(p!=j)/交换第p行与第j行for(k=0;k高斯列主元法解线性方程组for(n=j;n=0;i-)s=0.0;for(j=N-1;ji;j-)s=s+aij*xj;xi=(bi-s)/aii;cout文档视界经典四阶龙格库塔法解一阶微分方程组经典四阶龙格库塔法解一阶微分方程组高斯列主元法解线性方程组3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法讲明牛顿法主要思想是用xxF(x)F(

8、x)(k0,1,.).进行迭代。因而首先需要算出F(x)的雅可比矩阵F(x)，再求过F(x)求出它的逆F(x)，当它到达某个精度时即停止迭代。详细算法如下：首先设Pk已知：计算函数：计算雅可比矩阵：求线性方程组J(Pk)PF(Pk)的解P：计算下一点Pk1PkP重复上述经过F(Pk)f1(pk,qk)f2(pk,qk)3-1)J(Pk)f1(pk,qk)xf2(pk,qk)xf1(pk,qk)y3-2)3-3)(pk,qk)数值计算课程设计3.2牛顿法解非线性方程组流程图3.3牛顿法解非线性方程组程序代码#include#include#defineN2/#defineEpsilon0.000

9、1/#defineMax100/usingnamespacestd;constintN2=2*N;intmain()voidff(floatxxN,floatyyN);/非线性方程组中方程个数、未知量个数差向量1范数的上限最大迭代次数计算向量函数的因变量向量yyNvoidffjacobian(floatxxN,floatyyNN);/计算雅克比矩阵yyNNvoidinv_jacobian(floatyyNN,floatinvNN);/计算雅克比矩阵的逆矩阵invvoidnewdundiedai(floatx0N,floatinvNN,floatx1N);/由近似解向量x0计算近似解向量x1floaty0N,floatx0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errorno图3-1牛顿法解非线性方程组流程图