使用C语言解常微分方程 C.docx
《使用C语言解常微分方程 C.docx》由会员分享，可在线阅读，更多相关《使用C语言解常微分方程 C.docx（21页珍藏版）》请在冰豆网上搜索。
使用C语言解常微分方程C
解常微分方程
姓名：
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
0.
20
40
80
0.0005
0.000000214374
160
对比发现4阶精度远高于2阶精度
当我们细分到一定程度之后，我们发现误差的减小慢慢变小。
所以，若需要极高的精度的话会花费大量的计算资源。
II.解下面的方程组
我们选择M=1000来细分
运用3中代码，我们求解得
III.解下面高阶微分方程（震动方程）
运用3中代码，我们取步长h=,我们求解得
画出解y1,y2有，
IV.解洛伦兹方程
方程如下，使用不同的初始值，观察差别
分别使用初值
解查看附件
我们查看初始值和结束值。
t
x
y
z
x
y
z
0
5
5
5
5
5
我们发现尽管初值仅有很小的区别，但是终值有的巨大的差别（与相比）。
初值为
画出yz,xz,xy有，
初值为
画出yz,xz,xy有，
V.边值问题
我们考虑这样一个方程
使用3中代码求解有
详细答案参看附件。
我们看看几个均分点的值
t
y（t）打靶法
y（t）差分法
在算法上，打靶法计算量明显高于差分法，但是打靶法具有更高的普适性。
在进行，有解析解的方程求解时，发现在步长相同时，差分法具有更高的精度。
画出解的图有，
Shooting解法
有限差分解法
End.
代码：
文件
#include<>
#include<>
#include<>
#include""
/"
voidfree2DArray（double**a,intm）
{
inti;
for（i=0;ifree（a[i]）;
free（a）;
}
y[1];
y[1];
y[0]-4.*t*t-4.*t-2.;
}
voidfunv3（doublet,doubley[],doublefv[]）
{
fv[0]=y[1];
fv[1]=/*y[0]+/*y[2];
fv[2]=y[3];
fv[3]=/*y[2]+/*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）;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.};
=%,%\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.）=%,%\n",yv2[M-1],fabs（yv2[M-1]-sqrt（2.*exp（5.）+2.）））;
free（yv2）;
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（"","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,,11）;
fp=fopen（"","w+"）;
for（i=0;i<11;i++）
fprintf（fp,"%f\t%f\t%f\t%f\t%f\n",*i,yMa[i][0],yMa[i][1],yMa[i][2],yMa[i][3]）;
fclose（fp）;
free2DArray（yMa,11）;
yvLorenz,50.,M）;
fp=fopen（"","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]=;
yMa=RKSystem4th（funvLorenz,3,0.,yvLorenz,50.,M）;
fp=fopen（"","w+"）;
for（i=0;ifprintf（fp,"%f\t%f\t%f\n",yMa[i][0],yMa[i][1],yMa[i][2]）;
fclose（fp）;
tn=1.;
a=1.;
b=;
yMa=BVP_Shooting（funv4,t0,a,tn,b,M）;
fp=fopen（"","w+"）;
for（i=0;ifprintf（fp,"%f\t%f\n",t0+（tn-t0）/（M-1）*i,yMa[i][0]）;
fclose（fp）;
free2DArray（yMa,M）;
tn=1.;
a=1.;
b=;
yFD=BVP_FD（pfun,qfun,rfun,t0,a,tn,b,M）;
fp=fopen（"","w+"）;
for（i=0;ifprintf（fp,"%f\t%f\n",t0+（tn-t0）/（M-1）*i,yFD[i]）;
fclose（fp）;
free（yFD）;
return0;
}
文件
#include<>
#include""
#include<>
#include<>
/"
double*RungeKutta_Heum（doublefun（double,double）,doublet0,doubley0,doubletn,intM）
{
;
doublev[2]={a,-2.};
double**y;
doublex2,err;
doublef0,f1;
do{
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++;
}while（err>epsilon）;
returny;
}
/*double*BVP_FD（doublepfun（double）,doubleqfun（double）,doublerfun（double）,
doublet00,doublealpha,doubletn,doublebeta,intm）
{
-h*h*qfun（t）;
r[i]=h*h*rfun（t）;
t=t+h;
}
r[0]=r[0]-（h/2.*pfun（t0）+1.）*alpha;
-h/2.*pfun（tn））*beta;
t=t0;
for（i=0;i{
a[i]=1.+h/2.*pfun（t+h）;
c[i]=1.-h/2.*pfun（t）;
*（pfun（t0+h））-1.）*alpha;
r[m-1]=r[m-1]-（h/2*pfun（tn-h）-1）*beta;
d[0]=b[0];
u[0]=c[0]/b[0];
for（i=1;i{
d[i]=b[i]-a[i]*u[i-1];
if（d[i]==0）
{
return0;
}
u[i]=c[i]/d[i];
}
d[m-1]=b[m-1]-a[m-1]*u[m-2];
y[0]=r[0]/d[0];
for（i=1;i{
y[i]=（r[i]-a[i]*y[i-1]）/d[i];
}
x[m]=y[m-1];
for（i=m-2;i>=0;i--）
{
x[i+1]=y[i]-u[i]*x[i+2];
}
x[0]=alpha;
x[M-1]=beta;
returnx;
}
double*BVP_Lshooting（doublepfun（double）,doubleqfun（double）,doublerfun（double）,
doublet0,doublea,doubletn,doubleb,intM）
{
return0;
}
文件
#ifndefODE_HHHH
#defineODE_HHHH
//2ndorderRunge-Kuttamethodforsolvinginitialvalueproblem
double*RungeKutta_Heum（doublefun（double,double）,doublet0,doubley0,doubletn,intM）;
//4thorderRunge-Kuttamethodforsolvinginitialvalueproblem
double*RungeKutta4th（doublefun（double,double）,doublet0,doubley0,doubletn,intM）;
//ODEsystemsandhighorderODEsolver
double**RKSystem4th（voidfun（double,double[],double[]）,intN,doublet0,doubley0[],doubletn,intM）;
//ODEboundaryvalueproblemusingshootingmethod
double**BVP_Shooting（voidfun（double,double[],double[]）,doublet0,doublea,doubletn,doubleb,intM）;
//BVPusingfinitedifferencemethod
double*BVP_FD（doublepfun（double）,doubleqfun（double）,doublerfun（double）,
doublet0,doublea,doubletn,doubleb,intM）;
#endif
升级会员 