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

例题一
四阶龙格-库塔法的具体形式为：
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=Ae t – t -1, y=Be-t – t + 1
t=1时过（1,1）点的质点有
A=3/e , B=e
最后得此质点的迹线方程为：
x=3e t -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(derivative);
}
void main()
{
float t0; //范围上限
float t; //范围下限
float h; //步长
float nn; //计算出的点的个数,即迭代步数
int n; //对nn取整
float k1,k2,k3,k4,k11,k22,k33,k44;
float x[N],y[N]; //用于存放计算出的常微分方程数值解
float xx,yy; //精确值
float errorx,errory;//误差
int i=0,j;
//以下为函数的接口
printf("input t0:");
scanf("%f",&t0);
printf("input t:");
scanf("%f",&t);
printf("input x[0]:");
scanf("%f",&x[0]);
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++)
{
xx=3*exp(t0-1)-t0-1; //解析解表达式
yy=exp(1-t0)-t0+1;
errorx=x[i]-xx; //误差计算
errory=y[i]-yy;
printf("x[%f]=%f xx=%f errorx=%f\n",t0,x[i],xx,errorx);//结果输出printf("y[%f]=%f yy=%f errory=%f\n",t0,y[i],yy,errory);
k1=h*fx(t0,x[i]); //求K1
k2=h*fx((t0+h/2),(x[i]+k1*h/2)); //求K2
k3=h*fx((t0+h/2),(x[i]+k2*h/2)); //求K3
k4=h*fx((t0+h),(x[i]+k3*h)); //求K4
x[i+1]=x[i]+((k1+2*k2+2*k3+k4)/6); //求x[i+1]
k11=h*fy(t0,y[i]); //求K11
k22=h*fy((t0+h/2),(y[i]+k11*h/2)); //求K22
k33=h*fy((t0+h/2),(y[i]+k22*h/2)); //求K33
k44=h*fy((t0+h),(y[i]+k33*h)); //求K44
y[i+1]=y[i]+((k11+2*k22+2*k33+k44)/6); //求y[i+1]
t0+=float(h);
i++;
}
}
2.2 结果
input t0:1
input t:2
input x[0]:1
input y[0]:1
input h:0.2
nn=5.000000
n=5
ti xi xx errorx yi yy errory 1.0 1 1 0 1 1 0
1.2 1.428377 1.464208 -0.035831 0.588158 0.618731 -0.030573 1.4 1.984977
2.075475 -0.090498 0.217849 0.270320 -0.052471 1.6 2.695964 2.866357 -0.170393 -0.119071 -0.051189 -0.067882
1.8 3.592841 3.876624 -0.283783 -0.429147 -0.350671 -0.078476
2.0 4.713541 5.154847 -0.441306 -0.717643 -0.632121 -0.085522
图二：t=1至t=2质点的轨迹。