使用C语言解常微分方程 C.docx
- 文档编号:3964935
- 上传时间:2022-11-26
- 格式:DOCX
- 页数:21
- 大小:205.75KB
使用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;
打靶法
当err 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]); 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;i free(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;i fprintf(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;i fprintf(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;i fprintf(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;i fprintf(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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 使用C语言解常微分方程 使用 语言 微分方程