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

经典四阶龙格库塔法解一阶微分方程组 数值计算课程设计 PAGE \* MERGEFORMAT PAGE \* MERGEFORMAT 2 PAGE \* MERGEFORMAT PAGE \* MERGEFORMAT 1 1.经典四阶龙格库塔法解一阶微分方程组 1.1运用四阶龙格库塔法解一阶微分方程组算法分析 （1-1） （1-1） （1-2 （1-2） （1-3 （1-3） （1-4 （1-4） （1-5 （1-5） （1-6 （1-6） （1-7 （1-7） （1-8 （1-8） （1-9 （1-9） （1-10 （1-10） 经过循环计算由 推得 …… 每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。 1.2经典四阶龙格库塔法解一阶微分方程流程图 图1-1 经典四阶龙格库塔法解一阶微分方程流程图 1.3经典四阶龙格库塔法解一阶微分方程程序代码： #include iostream #include iomanip using namespace std; void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h) { double f1,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; } int main() { double f(double t,double x, double y); double g(double t,double x, double y); double initial[3],resu[3]; double a,b,H; double t,step; int i; 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; cout initial[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); } double f(double t,double x, double y) { double dx; dx=x+2*y; return(dx); } double g(double t,double x, double y) { double dy; dy=3*x+2*y; return(dy); } 1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示：