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

1、例题一 y'=t-y, 0t1 y0=0 四阶龙格-库塔法的具体形式为：K1=tn-yn K2=tn+t2-yn+K1×t2=0.9tn-yn +0.1 K3=tn+t2-yn+K2×t2=0.91tn-yn +0.09 K4=tn+t-yn+K3×t=0.818tn-yn +0.182yn+1=yn+t6K1+2K2+2K3+K4 1.1程序:/*e.g: y'=t-y,t0,1/* y(0)=0/*使用经典四阶龙格-库塔算法进行高精度求解/* y(i1)yi( K1+ 2*K2 +2*K3+ K4)/6/* K1=h*f(ti,yi)/* K2

2、=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 20float f(float d,float p) /要求解的微分方程的右部的函数 e.g: t-yfloat derivative; derivative=d-p;return(derivative);void main()float t0; /范围上限float t; /范围下限float h; /步长 float nn; /计算

3、出的点的个数,即迭代步数int n; /对nn取整 float k1,k2,k3,k4; float yN; /用于存放计算出的常微分方程数值解float yy; /精确值 float error;/误差int i=0,j;/以下为函数的接口printf("input t0:"); scanf("%f",&t0);printf("input t:"); scanf("%f",&t); printf("input y0:"); scanf("%f",&y

4、0);printf("input h:"); scanf("%f",&h);/ 以下为核心程序 nn=(t-t0)/h; printf("nn=%fn",nn);n=(int)nn;printf("n=%dn",n); for(j=0;j<=n;j+) yy=t0-1+exp(-t0); /解析解表达式 error=yi-yy; /误差计算 printf("y%f=%f yy=%f error=%fn",t0,yi,yy,error);/结果输出 k1=h*f(t0,yi); /求

5、K1 k2=h*f(t0+h/2),(yi+k1*h/2); /求K2 k3=h*f(t0+h/2),(yi+k2*h/2); /求K3 k4=h*f(t0+h),(yi+k3*h); /求K4 yi+1=yi+(k1+2*k2+2*k3+k4)/6); /求yi+1 t0+=float(h); i+; 1.2结果:input t0:0input t:1input y0:0input h:0.2nn=5.000000n=5tiyiyyerror00000.20.0197360.0187310.0010050.40.0748130.0703200.0044930.60.1583030.14881

6、20.0094910.80.2646350.2493290.0153061.00.3893310.3678790.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 + 1t=1时过（1,1）点的质点有 A=3/e , B=e最后得此质点的迹线方程为： x=3et -1 t -1, y=e-1-t t + 12.1 程序#include <stdio.h>#include <ma

7、th.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(derivative);void main()float t0; /范围上限float t; /范围下限float h; /步长 float nn; /计算出的点的个数,即迭代步数int n; /对nn取整 float k1,k2,k3,k4,k11,k22,k33,k44; float

8、xN,yN; /用于存放计算出的常微分方程数值解float xx,yy; /精确值 float errorx,errory;/误差int i=0,j;/以下为函数的接口printf("input t0:"); scanf("%f",&t0);printf("input t:"); scanf("%f",&t);printf("input x0:"); scanf("%f",&x0); printf("input y0:"); sca

9、nf("%f",&y0);printf("input h:"); scanf("%f",&h);/ 以下为核心程序 nn=(t-t0)/h; printf("nn=%fn",nn);n=int(nn);printf("n=%dn",n); for(j=0;j<=n;j+) xx=3*exp(t0-1)-t0-1; /解析解表达式 yy=exp(1-t0)-t0+1; errorx=xi-xx; /误差计算 errory=yi-yy; printf("x%f=%f

10、xx=%f errorx=%fn",t0,xi,xx,errorx);/结果输出 printf("y%f=%f yy=%f errory=%fn",t0,yi,yy,errory); k1=h*fx(t0,xi); /求K1 k2=h*fx(t0+h/2),(xi+k1*h/2); /求K2 k3=h*fx(t0+h/2),(xi+k2*h/2); /求K3 k4=h*fx(t0+h),(xi+k3*h); /求K4 xi+1=xi+(k1+2*k2+2*k3+k4)/6); /求xi+1 k11=h*fy(t0,yi); /求K11 k22=h*fy(t0+h/

11、2),(yi+k11*h/2); /求K22 k33=h*fy(t0+h/2),(yi+k22*h/2); /求K33 k44=h*fy(t0+h),(yi+k33*h); /求K44 yi+1=yi+(k11+2*k22+2*k33+k44)/6); /求yi+1 t0+=float(h); i+; 2.2 结果input t0:1input t:2input x0:1input y0:1input h:0.2nn=5.000000n=5tixixxerrorxyiyyerrory1.01101101.21.4283771.464208-0.0358310.5881580.618731-0.0305731.41.9849772.075475-0.0904980.2178490.270320-0.0524711.62