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

),,(1kkkyxtff?,)2,2,2(112ghyfhxhtffkkk????)2,2,2(223ghyfhxhtffkkk????),,(334hgyhfxhtffkkk????),,(1kkkyxtgg?)2,2,2(112ghyfhxhtggkkk????)2,2,2(223ghyfhxhtggkkk????),,(334hgyhfxhtggkkk????)22(6)22(64321143211gggghyyffffhxxkkkk????????????1k kt t h?? ?经过循环计算由推得……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为??NOh,一种折中方法是每次进行若干次函数求值,从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的,它适用于一般的应用,因为它非常精准,稳定,且易于编程。0 0 0, ,t x y????1 1 1 2 2 2, , , ,t x y t x y(1-1)(1-2)(1-3)(1-4)(1-5)(1-6)(1-7)(1-8)(1-9)(1-10)-:#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,数值计算课程设计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:";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;经典四阶龙格库塔法解一阶微分方程组4dx=x+2*y;return(dx);}doubleg(doublet,doublex,doubley){doubledy;dy=3*x+2*y;return(dy);}