使用C语言解常微分方程 C ODE.docx
《使用C语言解常微分方程 C ODE.docx》由会员分享，可在线阅读，更多相关《使用C语言解常微分方程 C ODE.docx（37页珍藏版）》请在冰点文库上搜索。
使用C语言解常微分方程CODE
解常微分方程
姓名：
Vincent
年级：
2010，学号：
1033****，组号：
5（小组），4（大组）
1.数值方法:
我们的实验目标是解常微分方程，其中包括几类问题。
一阶常微分初值问题，高阶常微分初值问题，常微分方程组初值问题，二阶常微分方程边值问题，二阶线性常微分方程边值问题。
对待上面的几类问题，我们分别使用不同的方法。
∙初值问题
使用龙格-库塔来处理
∙边值问题
用打靶法来处理
∙线性边值问题
有限差分法
初值问题
我们分别使用
∙二阶龙格-库塔方法
∙4阶龙格-库塔方法
来处理一阶常微分方程。
理论如下：
对于这样一个方程
当h很小时，我们用泰勒展开，
当我们选择正确的参数a[i][j],b[i][j]之后，就可以近似认为这就是泰勒展开。
对于二阶，我们有：
其中
经过前人的计算机经验，我们有，
选择A=1/2,B=1/2,则P=1,Q=1，于是又如下形式，这也使休恩方法的表达式。
所以我们称其为龙格（库塔）休恩方法
对于4阶龙格方法，我们有类似的想法，
我们使用前人经验的出的系数，有如下公式
对于高阶微分方程及微分方程组
我们用
∙4阶龙格-库塔方法来解
对于一个如下的微分方程组
我们可以认为是一个一阶向量微分方程，所以可以用龙格-库塔方法的向量形式解。
对于一个高阶的微分方程，形式如下:
我们可以构建出一个一阶的微分方程组，
令
则有
所以我们实际只要解一个微分方程组
其中初值为
使用4阶龙格-库塔方法，
使用这个向量形式的龙格-库塔方法我们便可就出方程的数值解。
边值问题
对于边值问题，我们分为两类
∙一般的边值问题
∙线性边值微分方程
一般的边值问题，我们是使用打靶法来求解，
对于这样一个方程
主要思路是，化成初值问题来求解。
我们已有
这样我们便可求出方程的解。
∙线性微分方程边值问题
对于这样的问题，我们可以使用一个更加高效的方法，有限差分法。
对于如下方程
我们对其进行差分
这样的话，我们的微分方程可以写成，
于是我们得到了个线性方程组
这样的话我们求解出x
对于上面的矩阵，我们可以分解出两个三角阵，然后回代求解便可。
我们用回代法可以直接求解
至此，我们便求出了目标方程的解
2.流程图
∙二阶龙格-库塔
对于i=0到M-1;
y[i+1]=y[i]+h/2*（fun（t,y[i]）+fun（t+h,y[i]+h*fun（t,y[i]）））;
returny;
∙4阶龙格-库塔
对于i=0到M-1;
y[i+1]=y[i]+h/6*（fun（t,y[i]）+2*fun（t+h/2,y[i]+h/2*fun（t,y[i]））+2*fun（t+h/2,y[i]+h/2*fun（t+h/2,y[i]+h/2*fun（t,y[i]）））+fun（t+h,y[i]+h*fun（t+h/2,y[i]+h/2*fun（t+h/2,y[i]+h/2*fun（t,y[i]）））））;
returny;
∙4阶龙格-库塔解方程组
对于k=0到M-1;
对于i=0到N;
fun（t,y[k],dy[0]）
对于i=0到N;
tempdy[j]=y[k][j]+h/2*dy[0][j];
fun（t+h/2,tempdy,dy[1]）;
对于i=0到N;
tempdy[j]=y[k][j]+h/2*dy[1][j];
fun（t+h/2,tempdy,dy[2]）;
对于i=0到N;
tempdy[j]=y[k][j]+h*dy[2][j];
fun（t+h,tempdy,dy[3]）;
y[k+1][i]=y[k][i]+h/6*（dy[0][i]+2*dy[1][i]+2*dy[2][i]+dy[3][i]）;
returny;
∙打靶法
当erry=RKSystem4th（fun,2,t0,u,tn,M）;
f0=y[M-1][0]-b;
u[1]=y[0][1];
y=RKSystem4th（fun,2,t0,v,tn,M）;
f1=y[M-1][0]-b;
v[1]=y[0][1];
x2=u[1]-f0*（v[1]-u[1]）/（f1-f0）;
u[1]=v[1];
v[1]=x2;
err=fabs（u[1]-v[1]）;
returny;
∙有限差分法
对i从0到m；求出，b,r,a,c
b[i]=2+h*h*qfun（t）;
r[i]=-h*h*rfun（t）;
a[i]=-h/2*pfun（t）-1;
c[i]=h/2*pfun（t）-1;
r[0]=r[0]-（（-h/2.）*（pfun（t0+h））-1.）*alpha;
r[m-1]=r[m-1]-（h/2*pfun（tn-h）-1）*beta;
求出d,u
d[0]=b[0];
u[0]=c[0]/b[0];
对i从1到m-1
d[i]=b[i]-a[i]*u[i-1];
u[i]=c[i]/d[i];
d[m-1]=b[m-1]-a[m-1]*u[m-2];
回代
y[0]=r[0]/d[0];
对于i从到m;
y[i]=（r[i]-a[i]*y[i-1]）/d[i];
x[m]=y[m-1];
对i从m–2到0
x[i+1]=y[i]-u[i]*x[i+2];
x[0]=alpha;
x[M-1]=beta;
returnx;
3.代码实现过程
查看附件
4.数值例子与分析
I.解下面的方程
求t=5时，y的值
使用3中代码执行得到
M
y（5）2阶方法解
y（5）4阶方法解
2阶方法误差
4阶方法误差
10
17.483960302367
17.287497436431
0.197366704867
0.000903838930
20
17.333302329975
17.286649150897
0.046708732474
0.000055553397
40
17.297973861450
17.286597041323
0.011380263949
0.000003443822
80
17.289403465806
17.286593811875
0.002809868305
0.000000214374
160
17.287291781639
17.286593610872
0.000698184138
0.000000013372
对比发现4阶精度远高于2阶精度
当我们细分到一定程度之后，我们发现误差的减小慢慢变小。
所以，若需要极高的精度的话会花费大量的计算资源。
II.解下面的方程组
我们选择M=1000来细分
运用3中代码，我们求解得
III.解下面高阶微分方程（震动方程）
运用3中代码，我们取步长h=0.02,我们求解得
画出解y1,y2有，
IV.解洛伦兹方程
方程如下，使用不同的初始值，观察差别
分别使用初值
解查看附件
我们查看初始值和结束值。
t
x
y
z
x
y
z
0
5
5
5
5.001
5
5
0.01
6.201596
0.657642
6.387359
6.202413
10.658526
6.387795
0.02
9.374014
17.452659
0.716395
9.374978
17.454007
10.717783
49.98
3.258276
3.369219
20.608133
-7.008305
-12.891724
11.712534
49.99
3.549458
4.588508
18.661441
-10.450006
-18.171700
16.666380
50.00
4.300485
6.279033
17.322649
-14.303620
-21.252383
26.374359
我们发现尽管初值仅有很小的区别，但是终值有的巨大的差别（与0.001相比）。
初值为
画出yz,xz,xy有，
初值为
画出yz,xz,xy有，
V.边值问题
我们考虑这样一个方程
使用3中代码求解有
详细答案参看附件。
我们看看几个均分点的值
t
y（t）打靶法
y（t）差分法
0.0
1.000000
1.000000
0.1
1.239224
1.239224
0.3
1.703017
1.703017
0.5
2.144261
2.144261
0.7
2.560903
2.560904
0.9
2.952845
2.952845
1.0
3.140000
3.140000
在算法上，打靶法计算量明显高于差分法，但是打靶法具有更高的普适性。
在进行，有解析解的方程求解时，发现在步长相同时，差分法具有更高的精度。
画出解的图有，
Shooting解法
有限差分解法
End.
代码：
文件main_ode.cpp
#include
#include
#include
#include"ode.h"
//#include"../array.h"
voidfree2DArray（double**a,intm）
{
inti;
for（i=0;ifree（a[i]）;
free（a）;
}
//y=f（t,y）,fun=f（t,y）
doublefun（doublet,doubley）
{
returnexp（t）/y;
}
voidfunv1（doublet,doubley[],doublefv[]）
{
//
//fv[0]=y[0]+2.*y[1];
//fv[1]=3*y[0]+2.*y[1];
//Lotka-Volterraequation
fv[0]=y[0]-y[0]*y[1]-0.1*y[0]*y[0];
fv[1]=y[0]*y[1]-y[1]-0.05*y[1]*y[1];
}
voidfunv2（doublet,doubley[],doublefv[]）
{
//y=x*x+x+1
fv[0]=y[1];
fv[1]=4.*y[0]-4.*t*t-4.*t-2.;
}
voidfunv3（doublet,doubley[],doublefv[]）
{
fv[0]=y[1];
fv[1]=-278.28/0.15*y[0]+110.57/0.15*y[2];
fv[2]=y[3];
fv[3]=-278.28/0.15*y[2]+110.57/0.15*y[0];
}
voidfunv4（doublet,doubley[],doublefv[]）
{
fv[0]=y[1];
fv[1]=-（t+1.）*y[1]+y[0]+（1.-t*t）*exp（-t）;
}
doublepfun（doublet）
{
return-（t+1.）;
}
doubleqfun（doublet）
{
return1.;
}
doublerfun（doublet）
{
return（1.-t*t）*exp（-t）;//-4.*t*t-4.*t-2.;
}
voidfunvLorenz（doublet,doubleyv[],doublefv[]）
{
doublex=yv[0],y=yv[1],z=yv[2];
fv[0]=10.*（-x+y）;
fv[1]=28.*x-y-x*z;
fv[2]=-8./3.*z+x*y;
}
intmain（intargc,char*argv[]）
{
inti,M;
doublea=0,b=0;
FILE*fp;
doublet0=0.,y0=2.,tn=5.,*yv,*yv2,**yMa,*yFD,yv0[2]={2.,1.},yvLorenz[3]={5.,5.,5.};
doubleyv3[4]={0.,1.,0.,1.};
//exactsolution:
y^2=2*exp（x）+2
//1storderODE
M=21;
yv=RungeKutta_Heum（fun,t0,y0,tn,M）;
printf（"y（5.）=%16.12f,%16.12f\n",yv[M-1],fabs（yv[M-1]-sqrt（2.*exp（5.）+2.）））;
free（yv）;
M=21;
yv2=RungeKutta4th（fun,t0,y0,tn,M）;
printf（"y（5.）=%16.12f,%16.12f\n",yv2[M-1],fabs（yv2[M-1]-sqrt（2.*exp（5.）+2.）））;
free（yv2）;
//ODEsystem
yMa=RKSystem4th（funv1,2,0.,yv0,30.,1000）;
/*yv0[0]=21.;
yv0[1]=-9.;
yMa=RKSystem4th（funv2,2,-5.,yv0,5.,3001）;*/
printf（"y1[30]=%f,y2[30]=%f\n",yMa[999][0],yMa[999][1]）;
/*printf（"erry1=%f,erry2=%f",4.*exp（4.）+2.*exp（-1.）-yMa[99][0],6.*exp（4.）-2.*exp（-1.）-yMa[99][1]）;
printf（"erry1=%f,erry2=%f",31-yMa[3000][0],11-yMa[3000][1]）;*/
free2DArray（yMa,100）;
yMa=RKSystem4th（funv1,2,0.,yv0,30.,1000）;
fp=fopen（"lv.dat","w+"）;
for（i=0;i<1000;i++）
fprintf（fp,"%f\t%f\n",yMa[i][0],yMa[i][1]）;
fclose（fp）;
free2DArray（yMa,1000）;
yMa=RKSystem4th（funv3,4,0.,yv3,0.2,11）;
fp=fopen（"fv3.dat","w+"）;
for（i=0;i<11;i++）
fprintf（fp,"%f\t%f\t%f\t%f\t%f\n",0.02*i,yMa[i][0],yMa[i][1],yMa[i][2],yMa[i][3]）;
fclose（fp）;
free2DArray（yMa,11）;
//Lorenzequ
M=1001;
yMa=RKSystem4th（funvLorenz,3,0.,yvLorenz,50.,M）;
fp=fopen（"Lorenz01.dat","w+"）;
for（i=0;ifprintf（fp,"%f\t%f\t%f\n",yMa[i][0],yMa[i][1],yMa[i][2]）;
fclose（fp）;
M=1001;
yvLorenz[0]=5.001;
yMa=RKSystem4th（funvLorenz,3,0.,yvLorenz,50.,M）;
fp=fopen（"Lorenz02.dat","w+"）;
for（i=0;ifprintf（fp,"%f\t%f\t%f\n",yMa[i][0],yMa[i][1],yMa[i][2]）;
fclose（fp）;
//highorderODE
//BVP
M=101;
t0=0.;
tn=1.;
a=1.;
b=3.14;
yMa=BVP_Shooting（funv4,t0,a,tn,b,M）;
fp=fopen（"Shoot.dat","w+"）;
for（i=0;ifprintf（fp,"%f\t%f\n",t0+（tn-t0）/（M-1）*i,yMa[i][0]）;
fclose（fp）;
free2DArray（yMa,M）;
//BVPFD
M=101;
t0=0.;
tn=1.;
a=1.;
b=3.14;
yFD=BVP_FD（pfun,qfun,rfun,t0,a,tn,b,M）;
fp=fopen（"yFD.dat","w+"）;
for（i=0;ifprintf（fp,"%f\t%f\n",t0+（tn-t0）/（M-1）*i,yFD[i]）;
fclose（fp）;
free（yFD）;
return0;
}
文件ode.cpp
#include
#include"ode.h"
#include
#include
//#include"../array.h"
double*RungeKutta_Heum（doublefun（double,double）,doublet0,doubley0,doubletn,intM）
{
//y（t+h）=y（t）+h/2*（f（t,y）+f（t+h,y+hf（t,y）））
doubleh=（tn-t0）/（M-1）;
doublet=t0;
//doubley[100]={0};
double*y;
y=（double*）malloc（M*sizeof（double））;
inti=0;
y[0]=y0;
for（i=0;i{
y[i+1]=y[i]+h/2*（fun（t,y[i]）+fun（t+h,y[i]+h*fun（t,y[i]）））;
t=t+h;
}
returny;
}
/*double*RungeKutta_EulerCauchy（doublefun（double,double）,doublea,doubleb,doubley0,intM）
{
}*/
double*RungeKutta4th（doublefun（double,double）,doublet0,doubley0,doubletn,intM）
{
doubleh=（tn-t0）/（M-1）;
doublet=t0;
//doubley[100]={0};
double*y;
y=（double*）malloc（M*sizeof（double））;
inti=0;
y[0]=y0;
for（i=0;i{
y[i+1]=y[i]+h/6*（fun（t,y[i]）+2*fun（t+h/2,y[i]+h/2*fun（t,y[i]））+
2*fun（t+h/2,y[i]+h/2*fun（t+h/2,y[i]+h/2*fun（t,y[i]）））+
fun（t+h,y[i]+h*fun（t+h/2,y[i]+h/2*fun（t+h/2,y[i]+h/2*fun（t,y[i]）））））;
t=t+h;
}
returny;
}
double**RKSystem4th（voidfun（double,double[],double[]）,intN,doublet0,doubley00[],doubletn,intM）
{
doubleh=（tn-t0）/（M-1）;
doublet=t0;
//doubley[100]={0};
double**y;
double**dy;
double*tempdy;
y=（double**）malloc（M*sizeof（double*））;
dy=（double**）malloc（4*sizeof（double*））;
tempdy=（double*）malloc（N*sizeof（double））;
inti=0,j=0,k=0;
for（i=0;i{
y[i]=（double*）malloc（N*sizeof（do
升级会员 