使用C语言解常微分方程 C ODEWord文档下载推荐.docx
《使用C语言解常微分方程 C ODEWord文档下载推荐.docx》由会员分享，可在线阅读，更多相关《使用C语言解常微分方程 C ODEWord文档下载推荐.docx（37页珍藏版）》请在冰点文库上搜索。
令
则有
所以我们实际只要解一个微分方程组
其中初值为
使用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阶龙格-库塔
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]）））））;
∙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]）;
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;
∙打靶法
当err<
epsilon
y=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]）;
∙有限差分法
对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
5
5.001
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有，
V.边值问题
我们考虑这样一个方程
使用3中代码求解有
详细答案参看附件。
我们看看几个均分点的值
y（t）打靶法
y（t）差分法
0.0
1.000000
0.1
1.239224
0.3
1.703017
0.5
2.144261
0.7
2.560903
2.560904
0.9
2.952845
1.0
3.140000
在算法上，打靶法计算量明显高于差分法，但是打靶法具有更高的普适性。
在进行，有解析解的方程求解时，发现在步长相同时，差分法具有更高的精度。
画出解的图有，
Shooting解法
有限差分解法
End.
代码：
文件main_ode.cpp
#include<
stdio.h>
stdlib.h>
math.h>
#include"
ode.h"
//#include"
../array.h"
voidfree2DArray（double**a,intm）
{
inti;
for（i=0;
i<
m;
i++）
free（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[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）;
yv2=RungeKutta4th（fun,t0,y0,tn,M）;
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）;
*/
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]）;
31-yMa[3000][0],11-yMa[3000][1]）;
free2DArray（yMa,100）;
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"
"
for（i=0;
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]）;
free2DArray（yMa,11）;
//Lorenzequ
M=1001;
yMa=RKSystem4th（funvLorenz,3,0.,yvLorenz,50.,M）;
Lorenz01.dat"
M;
%f\t%f\t%f\n"
yMa[i][0],yMa[i][1],yMa[i][2]）;
yvLorenz[0]=5.001;
yMa=RKSystem4th（funvLorenz,3,0.,yvLorenz,50.,M）;
Lorenz02.dat"
yMa[i][0],yMa[i][1],yMa[i][2]）;
//highorderODE
//BVP
M=101;
t0=0.;
tn=1.;
a=1.;
b=3.14;
yMa=BVP_Shooting（funv4,t0,a,tn,b,M）;
Shoot.dat"
t0+（tn-t0）/（M-1）*i,yMa[i][0]）;
free2DArray（yMa,M）;
//BVPFD
yFD=BVP_FD（pfun,qfun,rfun,t0,a,tn,b,M）;
yFD.dat"
t0+（tn-t0）/（M-1）*i,yFD[i]）;
free（yFD）;
return0;
文件ode.cpp
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;
i<
M-1;
{
y[i+1]=y[i]+h/2*（fun（t,y[i]）+fun（t+h,y[i]+h*fun（t,y[i]）））;
t=t+h;
}
/*double*RungeKutta_EulerCauchy（doublefun（double,double）,doublea,doubleb,doubley0,intM）
}*/
double*RungeKutta4th（doublefun（double,double）,doublet0,doubley0,doubletn,intM）
doubleh=（tn-t0）/（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]）））））;
double**RKSystem4th（voidfun（double,double[],double[]）,intN,doublet0,doubley00[],doubletn,intM）
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;
M;
y[i]=（double*）malloc（N*sizeof（do
升级会员 