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