四阶龙格库塔法的编程(赵).docVIP

例题一 四阶龙格-库塔法的具体形式为： 1.1程序: /*e.g: y=t-y,t∈[0,1] /*　　 y(0)=0 /*使用经典四阶龙格-库塔算法进行高精度求解　 /* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6 /* K1=h*f(ti,yi) /* K2=h*f(ti+h/2,yi+K1*h/2) /* K3=h*f(ti+h/2,yi+K2*h/2) /* K4=h*f(ti+h,yi+K3*h) */ #include stdio.h #include math.h #define N 20 float f(float d,float p) //要求解的微分方程的右部的函数 e.g: t-y { float derivative; derivative=d-p; return(derivative); } void main() { float t0; //范围上限 float t; //范围下限 float h; //步长 float nn; //计算出的点的个数,即迭代步数 int n; //对nn取整 float k1,k2,k3,k4; float y[N]; //用于存放计算出的常微分方程数值解 float yy; //精确值 float error;//误差 int i=0,j; //以下为函数的接口 printf(input t0:); scanf(%f,t0); printf(input t:); scanf(%f,t); printf(input y[0]:); scanf(%f,y[0]); printf(input h:); scanf(%f,h); // 以下为核心程序 nn=(t-t0)/h; printf(nn=%f\n,nn); n=(int)nn; printf(n=%d\n,n); for(j=0;j=n;j++) { yy=t0-1+exp(-t0); //解析解表达式 error=y[i]-yy; //误差计算 printf(y[%f]=%f yy=%f error=%f\n,t0,y[i],yy,error);//结果输出 k1=h*f(t0,y[i]); //求K1 k2=h*f((t0+h/2),(y[i]+k1*h/2)); //求K2 k3=h*f((t0+h/2),(y[i]+k2*h/2)); //求K3 k4=h*f((t0+h),(y[i]+k3*h)); //求K4 y[i+1]=y[i]+((k1+2*k2+2*k3+k4)/6); //求y[i+1] t0+=float(h); i++; } } 1.2结果: input t0:0 input t:1 input y[0]:0 input h:0.2 nn=5.000000 n=5 ti yi yy error 0 0 0 0 0.2 0.019736 0.018731 0.001005 0.4 0.074813 0.070320 0.004493 0.6 0.158303 0.148812 0.009491 0.8 0.264635 0.249329 0.015306 1.0 0.389331 0.367879 0.021452 图一 例题二(见《高等流体力学》P45页 例1.6) 给定速度场u=x+t, v=-y-t, 求t=1时过（1,1）点的质点的迹线。 解答：由迹线微分方程dx/dt= x+t, dy/dt=-y-t 积分得 x=Aet – t -1, y=Be-t – t + 1 t=1时过（1,1）点的质点有 A=3/e , B=e 最后得此质点的迹线方程为： x=3et -1– t -1, y=e-1-t – t + 1 2.1 程序 #include stdio.h #include math.h #define N 20 //定义微分方程 float fx(float dd,float pp) { float der; der=dd+pp; return(der); } float fy(float d,float p) { float derivative; derivative=-d-p; return(derivati