c语言 常微分方程 库,使用C语言解常微分方程CODE.docx

使用C语言解常微分方程CODE.docx

.解常微分方程姓名 Vincent年级 2010，学号 1033 ，组号 5(小组)，4(大组)1. 数值方法 我们的实验目标是解常微分方程，其中包括几类问题。一阶常微分初值问题，高阶常微分初值问题，常微分方程组初值问题，二阶常微分方程边值问题，二阶线性常微分方程边值问题。对待上面的几类问题，我们分别使用不同的方法。 初值问题使用 龙格 -库塔 来处理 边值问题用打靶法来处理 线性边值问题有限差分法初值问题我们分别使用 二阶 龙格 - 库塔 方法 4 阶 龙格 - 库塔 方法来处理一阶常微分方程。理论如下对于这样一个方程y tf t, y当 h 很小时，我们用泰勒展开，k1hf tk , yk k2hf tka1h, ykLkihf tkai h, yk当我们选择正确的参数对于二阶，我们有b1 1k1i1hbij k j j1aij,bij 之后，就可以近似认为这就是泰勒展开。y thy th Af 0Bf 1其中f0f t, yf1f tPh, yQhf0 经过前人的计算机经验，我们有，选择 A1/2,B1/2,则 P1,Q1，于是又如下形式，这也使休恩方法的表达式。所以我们称其为龙格(库塔)休恩方法hythytf t , yf th, yhf t, y2对于 4 阶龙格方法，我们有类似的想法，我们使用前人经验的出的系数，有如下公式yn 1ynh2k2 2k3 k4 ,k16k1f tn , yn ,hhk2f tn2 , yn2 k1 ,k3f tnh , ynh k2 ,22k4f tnh, ynhk3 .对于高阶微分方程及微分方程组我们用 4 阶龙格 -库塔方法来解对于一个如下的微分方程组y1f1 t, y1,L , yn ,Lynfn t , y1,L , yn .y1 t0 y1,0 ,Lyn t0 yn ,0我们可以认为是一个一阶向量微分方程，所以可以用龙格-库塔方法的向量形式解。对于一个高阶的微分方程，形式如下y n xfL n 1,t, y, y , , yyt0 0 ,Ly n1 t0 n 1我们可以构建出一个一阶的微分方程组，令y1 ty t ,Lyn 1tyn1 t则有yn1 xf t, y, y1 ,L , yn 1 ,yn2 xyn 1 x,My2 xy2 x,y xy1 xyt0 0 ,其中，初值为y1 t0 1 ,Lyn 1t0 n1所以我们实际只要解一个微分方程组y1f1 t, y1,L, yn ,y1t0 y1,0 ,L其中初值为Lynfn t , y1,L, yn .yn t0 yn,0使用 4 阶龙格 -库塔方法，对于 k1,n有,ym 1ymhfk,12 fk ,22 fk ,3fk,4kk6其中 , h 是步长， m是递归的点。fk ,1fk tm, y1m, y2m ynm ffth , ymh f, ymh f ymh fk ,2km2121,1222,1n2n,1fk ,3fktmh , y1mhf1,2 , y2mhf2,2 ynmhfn,2 2222fk ,2fk tmh, y1mhf1,3 , y2mhf2,3 ynmhfn,3 使用这个向量形式的龙格-库塔方法我们便可就出方程的数值解。边值问题对于边值问题，我们分为两类 一般的边值问题 线性边值微分方程一般的边值问题，我们是使用打靶法来求解，对于这样一个方程x f t , x, x, atb边界为， xa, xb.主要思路是，化成初值问题来求解。我们已有xa,我们估计 xamk利用初值问题方法求解出xbsk我们用割线法迭代，使得在进行一定数量的步骤后，sk ;这样我们便可求出方程的解。 线性微分方程边值问题对于这样的问题，我们可以使用一个更加高效的方法，有限差分法。对于如下方程x t pt xtqt xt r t ta, bxa, xb其中 p, q, r 是 t 的函数我们对其进行差分tnb aa nh, hNh 是步长当 h 足够小时，我们有x tn xn 12xnxn 1h2x tn xn 1xn 1 .2h这样的话，我们的微分方程可以写成，x2xxh2pt xn 1xn 1qtxnr t0n 1nn 1n2hnn1h2qnxn1h2L2 pn xn 12 h2 pn xn 1hrn , n 1, , N 1x0, xN于是我们得到了个线性方程组.2hr2h q112p10h2hx112 p22 hq212 p2x2h pnh pn12 h2 qn1xn22xN 2h pN 22 h2qN 21 h pN 21xN 1r22h pN 12 h2qN 1012其中hhe0 12 p1 , en 12 pN 1这样的话我们求解出x对于上面的矩阵，我们可以分解出两个三角阵，然后回代求解便可。b1c1x1f1a2b2c2x2f2OOOMMA x fan 1bn 1cn 1xn 1f n 1anbnxnf n其中 A 可以写成d11u1l 2d21u2AOOOln 1dn 11un 1l n dn1其中1l kak , k2,L,n2d1b1, u1c1 / b1 ,3dkbkak uk 1for k2,L, n1ukck / dk4 dn bn anun 1我们用回代法可以直接求解.2h r1e02h r22h rnh2rN 2h2 rN 1enAxfLUxfLy f ,Ux y回代求出 yy1f1 / d1yk fkak yk 1 / dk , k 2,L , n再次回代，求出xxnynxkykLuk yk 1, k n 1, ,1至此，我们便求出了目标方程的解2. 流程图 二阶龙格 - 库塔对于 i 0 到 M-1;yi1 yi h / 2 * funt, yi funt h, yi h*funt, yi;return y; 4 阶龙格 -库塔对于 i 0 到 M-1;yi 1 yi h / 6 * funt, yi 2 * funt h / 2, yi h / 2 * funt, yi 2 * funt h / 2, yi h / 2 * funt h / 2, yi h / 2 * funt, yi funt h, yih*funt h / 2, yi h / 2 * funt h / 2, yi h / 2 *funt, yi;return y; 4 阶龙格 -库塔解方程组对于 k 0 到 M-1;对于 i 0 到 N;funt, yk, dy0对于 i 0 到 N;tempdyj ykj h / 2 * dy0j;funt h / 2, tempdy, dy1;对于 i 0 到 N;tempdyj ykj h / 2 * dy1j;funt h / 2, tempdy, dy2;对于 i 0 到 N;tempdyj ykj h * dy2j;funt h, tempdy, dy3;yk1i yki h / 6 * dy0i 2 * dy1i 2 * dy2i dy3i;return y; 打靶法当 errepsilony RKSystem4th fun, 2, t0, u, tn, M;f0 yM-10 - b;u1 y01;y RKSystem4th fun, 2, t0, v, tn, M;f1 yM-10 - b;v1 y01;x2 u1 - f0 * v1-u1/f1-f0;u1 v1;v1 x2;err fabsu1 - v1;return y; 有限差分法对 i 从 0 到 m；求出， b,r,a,cbi 2 h*h*qfunt; ri -h*h*rfunt;ai -h / 2 * pfunt - 1;ci h / 2 * pfunt - 1;r0 r0 - -h / 2. * pfunt0h - 1.*alpha;rm - 1 rm - 1 - h / 2 * pfuntn - h - 1*beta;求出 d,ud0 b0;u0 c0 / b0;对 i 从 1 到 m - 1di bi - ai * ui - 1;ui ci / di;dm - 1 bm - 1 - am - 1 * um - 2;回代y0 r0 / d0;对于 i 从 到 m;yi ri - ai * yi - 1 / di;xm ym -1;对 i 从 m 2 到 0 xi 1 yi - ui * xi 2;x0 alpha;xM - 1 beta;return x;3. 代码实现过程查看附件4. 数值例子与分析I. 解下面的方程y et / y,y02.求 t5 时， y 的值使用 3 中代码执行得到My5 2 阶方法解y54 阶方法解2 阶方法误差4 阶方法误差1017.48396030236717.2874974364310.1973667048670.0009038389302017.33330232997517.2866491508970.0467087324740.0000555533974017.29797386145017.2865970413230.0113802639490.0000034438228017.28940346580617.2865938118750. 0028098683050. 00000021437416017.28729178163917.2865936108720.0006981841380.000000013372对比发现4 阶精度远高于2 阶精度当我们细分到一定程度之后，我们发现误差的减小慢慢变小。所以，若需要极高的精度的话会花费大量的计算资源。II. 解下面的方程组x f t, x, yxxy101 x2 ,y gt, x, yxyy201 y2 ,x02, y01,t0,30.我们选择M1000 来细分运用 3 中代码，我们求解得x 300.990034y 300.842214III. 解下面高阶微分方程 震动方程 0.15 y1278.78y1110.57y20,0.15 y2278.78y2110.57y10,y1 0y2 00,y1 0y2 01,t0,0.2.运用 3 中代码，我们取步长h0.02, 我们求解得txtxtytyt0.0000000.0000001.0000000.0000001.0000000.0200000.0185090.7847200.0185090.7847200.0400000.0290490.2327450.0290490.2327450.0600000.027103-0.4185200.027103-0.4185200.0800000.013522-0.8893150.013522-0.8893150.100000-0.005849-0.977698-0.005849-0.9776980.120000-0.022687-0.646168-0.022687-0.6461680.140000-0.029763-0.037571-0.029763-0.0375710.160000-0.0240510.586445-0.0240510.586445画出解 y1,y2 有，IV. 解洛伦兹方程方程如下，使用不同的初始值，观察差别dx10 x10 y,dtdy28 xy xz,dtdz8xy.dt3 z分别使用初值 x0 , y0 , z0 5,5,5, x0 , y0 , z0 5.001,5,5解查看附件我们查看初始值和结束值。 x0 , y0 , z0 5,5,5 x0 , y0 , z0 5.001,5,5txyzxyz05555.001550.016.2015960.6576426.3873596.20241310.6585266.3877950.029.37401417.4526590.7163959.37497817.45400710.71778349.983.2582763.36921920.608133-7.008305-12.89172411.71253449.993.5494584.58850818.66144-10.450006-18.17170016.666380150.004.3004856.27903317.322649-14.303620-21.25238326.374359我们发现尽管初值仅有很小的区别，但是终值有的巨大的差别(与0.001 相比)。初值为 x0 , y0 , z0 5,5,5画出 yz,xz,xy 有，初值为 x0 , y0 , z0 5.001,5,5画出 yz,xz,xy 有，V. 边值问题我们考虑这样一个方程y t1y 2y1t 2 e t , 0t1y01; y13.14使用 3 中代码求解有详细答案参看附件。我们看看几个均分点的值tyt 打靶法yt 差分法0.01.0000001.0000000.11.2392241.2392240.31.7030171.7030170.52.1442612.1442610.72.5609032.5609040.92.9528452.9528451.03.1400003.140000在算法上，打靶法计算量明显高于差分法，但是打靶法具有更高的普适性。在进行，有解析解的方程求解时，发现在步长相同时，差分法具有更高的精度。画出解的图有，Shooting 解法有限差分解法End.代码文件 main_ode.cppinclude stdio.hinclude stdlib.hinclude math.hinclude ode.hinclude /array.hvoid free2DArraydouble a, int mint i;for i0; im; i freeai;freea; y ft,y, fun ft,ydouble fundouble t,double yreturn expt/y;void funv1double t,double y,double fv fv0 y0 2.*y1; fv13*y0 2.*y1; Lotka-Volterra equationfv0 y0 - y0*y1 - 0.1 *y0*y0;fv1 y0*y1 - y1 - 0.05*y1*y1;void funv2double t,double y,double fv yx*xx1 fv0 y1;fv1 4.*y0 - 4.*t*t - 4.*t - 2.;void funv3double t, double y, double fvfv0 y1;fv1 -278.28 / 0.15*y0 110.57 / 0.15*y2;fv2 y3;fv3 -278.28 / 0.15*y2 110.57 / 0.15*y0;void funv4double t, double y, double fvfv0 y1;fv1 -t 1.*y1 y0 1. - t*t*exp-t;double pfundouble treturn -t1.;double qfundouble treturn 1.;double rfundouble treturn 1. - t*t*exp-t; -4.*t*t - 4.*t - 2.;void funvLorenzdouble t,double yv,double fvdouble xyv0, yyv1, zyv2;fv0 10.*-xy;fv1 28.*x - y - x*z;fv2 -8./ 3.*z x*y;int main int argc, char *argv int i,M;double a 0, b 0;FILE *fp;double t00.,y02., tn5., *yv,*yv2, yMa, *yFD, yv022.,1., yvLorenz35.,5.,5.; double yv34 0., 1., 0., 1. ; exact solution y2 2*expx2 1st order ODEM21;yv RungeKutta_Heum fun, t0, y0, tn, M;printfy5.16.12f, 16.12f n, yvM-1, fabsyvM-1-sqrt2.*exp5.2.; freeyv;M21;yv2 RungeKutta4th fun, t0, y0, tn, M;printfy5.16.12f, 16.12f n, yv2M-1, fabsyv2M-1-sqrt2.*exp5.2.;freeyv2; ODE systemyMa RKSystem4th funv1, 2, 0., yv0, 30., 1000;/*yv00 21.;yv01 -9.;yMa RKSystem4th funv2, 2, -5., yv0, 5., 3001;*/printf y130f, y230f n, yMa9990,yMa9991;/*printferry1f,erry2f, 4. * exp4. 2. * exp-1. - yMa990, 6. * exp4. - 2.*exp-1. -yMa991;printferry1f,erry2f, 31 - yMa30000, 11 - yMa30001;*/ free2DArrayyMa,100;yMa RKSystem4th funv1, 2, 0., yv0, 30., 1000;fplv.dat,w;fori0;i1000;ifprintffp,ft fn,yMai0,yMai1;fclosefp;free2DArrayyMa,1000;yMa RKSystem4thfunv3, 4, 0., yv3, 0.2, 11;fp fv3.dat, w;for i 0; i11; ifprintffp,ftftftftfn,0.02*i,yMai0,yMai1,yMai2,yMai3;fclosefp;free2DArrayyMa, 11; Lorenz equM