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

(1-1),(1-2)(1-3)(1-4)(1-5)(1-6)(1-7)(1-8)(1-9)(1-10)经过循环计算由推得……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为,一种折中方法是每次进行若干次函数求值,从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的,它适用于一般的应用,因为它非常精准,稳定,且易于编程。-:#include<iostream>#include<iomanip>usingnamespacestd;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,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;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:";cin>>initial[0]>>initial[1]>>initial[2];cout<<"输入所求微分方程组的微分区间[a,b]:";cin>>a>>b;cout<<"输入所求微分方程组所分解子区间的个数step:";cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision(10);H=(b-a)/step;cout<<initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;for(i=0;i<step;i++){RK4(f,g,initial,resu,H);cout<<resu[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;dx=x+2*y;return(dx);}doubleg(doublet,doublex,doubley){doubledy;dy=3*x+2*y;return(dy);}:应用所编写程序计算所给例题:其中初值为 求解区间为[0,]。计算结果为:图1-:fork=1ton-1do fori=k+1ton l<=a(i,k)/a(k,k) forj=ktondo a(i,j)<=a(i,j)-l*a(k,j)end%endofforj b(i)<=b(i)-l*b(k)end%endoffori end%endoffork最后得到A,b可以构成上三角线性方程组接着使用回代法求解上三角线性方程组因为高斯消元要求a(k,k)≠0(k=1,2,3……n-1)这就需要对高斯消元过程进行完善,即使用高斯列主元法:其步骤为:①找主元:计算,并记录其所在行r,②交换第r行与第k行;③以第k行为工具行处理以下各行,使得从第k列的第k+