使用C語(yǔ)言解常微分方程CODE

1、.解常微分方程姓名： Vincent年級(jí)： 2010，學(xué)號(hào)： 1033* ，組號(hào)： 5（小組），4（大組）1. 數(shù)值方法 :我們的實(shí)驗(yàn)?zāi)繕?biāo)是解常微分方程，其中包括幾類問題。一階常微分初值問題，高階常微分初值問題，常微分方程組初值問題，二階常微分方程邊值問題，二階線性常微分方程邊值問題。對(duì)待上面的幾類問題，我們分別使用不同的方法。? 初值問題使用 龍格 -庫(kù)塔 來處理? 邊值問題用打靶法來處理? 線性邊值問題有限差分法初值問題我們分別使用? 二階 龍格 - 庫(kù)塔 方法? 4 階 龍格 - 庫(kù)塔 方法來處理一階常微分方程。理論如下：對(duì)于這樣一個(gè)方程y (t)f (t, y)當(dāng) h 很小時(shí)，我們用泰

2、勒展開，k1hf (tk , yk )k2hf (tka1h, ykLkihf (tkai h, yk當(dāng)我們選擇正確的參數(shù)對(duì)于二階，我們有：b1 1k1)i1hbij k j )j1aij,bij 之后，就可以近似認(rèn)為這就是泰勒展開。y( th)y( t)h Af 0Bf 1其中f0f (t, y)f1f (tPh, yQhf0 )經(jīng)過前人的計(jì)算機(jī)經(jīng)驗(yàn)，我們有，.選擇 A=1/2,B=1/2,則 P=1,Q=1，于是又如下形式，這也使休恩方法的表達(dá)式。所以我們稱其為龍格（庫(kù)塔）休恩方法hy(th)y(t)f (t , y)f (th, yhf (t, y)2對(duì)于 4 階龍格方法，我們有類似的想

3、法，我們使用前人經(jīng)驗(yàn)的出的系數(shù)，有如下公式y(tǒng)n 1ynh2k2 2k3 k4 ),(k16k1f ( tn , yn ),hhk2f (tn2 , yn2 k1 ),k3f (tnh , ynh k2 ),22k4f (tnh, ynhk3 ).對(duì)于高階微分方程及微分方程組我們用? 4 階龍格 -庫(kù)塔方法來解對(duì)于一個(gè)如下的微分方程組y1f1 (t, y1,L , yn ),Lynfn (t , y1,L , yn ).y1 (t0 )y1,0 ,Lyn (t0 ) yn ,0我們可以認(rèn)為是一個(gè)一階向量微分方程，所以可以用龍格-庫(kù)塔方法的向量形式解。對(duì)于一個(gè)高階的微分方程，形式如下:y( n )

4、( x)fL( n 1),t, y, y , , yy(t0 )0 ,Ly( n1) (t0 )n 1我們可以構(gòu)建出一個(gè)一階的微分方程組，令y1 (t)y (t ),Lyn 1(t)y(n1) (t)則有.yn1( x)f t, y, y1 ,L , yn 1 ,yn2 ( x)yn 1 ( x),My2 ( x)y2 (x),y (x)y1( x)y(t0 )0 ,其中，初值為y1 (t0 )1 ,Lyn 1(t0 )n1所以我們實(shí)際只要解一個(gè)微分方程組y1f1 (t, y1,L, yn ),y1(t0 )y1,0 ,L其中初值為L(zhǎng)ynfn (t , y1,L, yn ).yn (t0 )y

5、n,0使用 4 階龍格 -庫(kù)塔方法，對(duì)于 k1,n有,ym 1ymhfk,12 fk ,22 fk ,3+fk,4kk6其中 , h 是步長(zhǎng)， m是遞歸的點(diǎn)。fk ,1fk (tm, y1m, y2m,., ynm )ff(th , ymh f, ymh f,., ymh f)k ,2km2121,1222,1n2n,1fk ,3fk(tmh , y1mhf1,2 , y2mhf2,2 ,., ynmhfn,2 )2222fk ,2fk (tmh, y1mhf1,3 , y2mhf2,3 ,., ynmhfn,3 )使用這個(gè)向量形式的龍格-庫(kù)塔方法我們便可就出方程的數(shù)值解。邊值問題對(duì)于邊值問題

6、，我們分為兩類? 一般的邊值問題? 線性邊值微分方程一般的邊值問題，我們是使用打靶法來求解，對(duì)于這樣一個(gè)方程x f (t , x, x), atb邊界為， x(a), x(b).主要思路是，化成初值問題來求解。我們已有.x(a),我們估計(jì) x(a)mk利用初值問題方法求解出x(b)sk我們用割線法迭代，使得在進(jìn)行一定數(shù)量的步驟后，sk ;這樣我們便可求出方程的解。? 線性微分方程邊值問題對(duì)于這樣的問題，我們可以使用一個(gè)更加高效的方法，有限差分法。對(duì)于如下方程x (t )p(t) x(t)q(t )x(t )r (t )ta, bx(a), x(b)其中 p, q, r 是 t 的函數(shù)我們對(duì)其進(jìn)

7、行差分tnb aa nh, hNh 是步長(zhǎng)當(dāng) h 足夠小時(shí)，我們有x (tn )xn 12xnxn 1h2x (tn )xn 1xn 1 .2h這樣的話，我們的微分方程可以寫成，x2xxh2p(t) xn 1xn 1q(t)xnr (t)0n 1nn 1n2hnn1h2qnxn1h2L2 pn xn 12 h2 pn xn 1hrn , n 1, , N 1x0, xN于是我們得到了個(gè)線性方程組.2hr2h q112p10h2hx112 p22 hq212 p2x2h pnh pn12 h2 qn1xn22xN 2h pN 22 h2qN 21 h pN 21xN 1r22h pN 12 h2

8、qN 1012其中hhe0( 12 p1 ) , en( 12 pN 1)這樣的話我們求解出x對(duì)于上面的矩陣，我們可以分解出兩個(gè)三角陣，然后回代求解便可。b1c1x1f1a2b2c2x2f2OOOMMA x fan 1bn 1cn 1xn 1f n 1anbnxnf n其中 A 可以寫成d11u1l 2d21u2AOOOln 1dn 11un 1l n dn1其中(1)l kak , k2,L,n(2)d1b1, u1c1 / b1 ,(3)dkbkak uk 1for k2,L, n1ukck / dk(4) dn bn anun 1我們用回代法可以直接求解.2h r1e02h r22h r

9、nh2rN 2h2 rN 1en.AxfLUxfLy f ,Ux y回代求出 yy1f1 / d1yk( fkak yk 1) / dk , k 2,L , n再次回代，求出xxnynxkykLuk yk 1, k n 1, ,1至此，我們便求出了目標(biāo)方程的解2. 流程圖? 二階龍格 - 庫(kù)塔對(duì)于 i = 0 到 M-1;yi+1 = yi + h / 2 * (fun(t, yi) + fun(t + h, yi + h*fun(t, yi);return y;? 4 階龍格 -庫(kù)塔對(duì)于 i = 0 到 M-1;yi + 1 = yi + h / 6 * (fun(t, yi) + 2 *

10、fun(t + h / 2, yi + h / 2 * fun(t, yi) + 2 * fun(t + h / 2, yi + h / 2 * fun(t + h / 2, yi + h / 2 * fun(t, yi) +fun(t + h, yi+h*fun(t + h / 2, yi + h / 2 * fun(t + h / 2, yi + h / 2 *fun(t, yi);return y;? 4 階龍格 -庫(kù)塔解方程組對(duì)于 k= 0 到 M-1;對(duì)于 i= 0 到 N;fun(t, yk, dy0)對(duì)于 i= 0 到 N;tempdyj = ykj + h / 2 * dy0j

11、;fun(t + h / 2, tempdy, dy1);對(duì)于 i= 0 到 N;tempdyj = ykj + h / 2 * dy1j;fun(t + h / 2, tempdy, dy2);對(duì)于 i= 0 到 N;tempdyj = ykj + h * dy2j;fun(t + h, tempdy, dy3);yk+1i = yki + h / 6 * (dy0i + 2 * dy1i + 2 * dy2i + dy3i);return y;? 打靶法當(dāng) errepsilony = RKSystem4th( fun, 2, t0, u, tn, M);f0 = yM-10 - b;u1

12、= 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 = fabs(u1 - v1);return y;? 有限差分法對(duì) i 從 0 到 m；求出， b,r,a,cbi = 2 + h*h*qfun(t); ri = -h*h*rfun(t);ai = -h / 2 * pfun(t) - 1;ci = h / 2 * pfun(t) - 1;r0 = r0 - (-h / 2. )*( pfun(t0+h)

13、- 1.)*alpha;rm - 1 = rm - 1 - (h / 2 * pfun(tn - h) - 1)*beta;求出 d,ud0 = b0;u0 = c0 / b0;對(duì) i 從 1 到 m - 1di = bi - ai * ui - 1;ui = ci / di;dm - 1 = bm - 1 - am - 1 * um - 2;回代y0 = r0 / d0;對(duì)于 i 從 到 m;yi = (ri - ai * yi - 1) / di;xm = ym -1;對(duì) i 從 m 2 到 0xi + 1 = yi - ui * xi + 2;x0 = alpha;xM - 1 = be

14、ta;return x;3. 代碼實(shí)現(xiàn)過程查看附件4. 數(shù)值例子與分析I. 解下面的方程y et / y,y(0)2.求 t=5 時(shí)， y 的值使用 3 中代碼執(zhí)行得到.My(5) 2 階方法解y(5)4 階方法解2 階方法誤差4 階方法誤差1017.48396030236717.2874974364310.1973667048670.0009038389302017.33330232997517.2866491508970.0467087324740.0000555533974017.29797386145017.2865970413230.0113802639490.00000344382

15、28017.28940346580617.2865938118750. 0028098683050. 00000021437416017.28729178163917.2865936108720.0006981841380.000000013372對(duì)比發(fā)現(xiàn)4 階精度遠(yuǎn)高于2 階精度當(dāng)我們細(xì)分到一定程度之后，我們發(fā)現(xiàn)誤差的減小慢慢變小。所以，若需要極高的精度的話會(huì)花費(fèi)大量的計(jì)算資源。II. 解下面的方程組x f (t, x, y)xxy101 x2 ,y g(t, x, y)xyy201 y2 ,x(0)2, y(0)1,t0,30.我們選擇M=1000 來細(xì)分運(yùn)用 3 中代碼，我們求解得x 3

16、00.990034y 300.842214III. 解下面高階微分方程 (震動(dòng)方程 )0.15 y1278.78y1110.57y20,0.15 y2278.78y2110.57y10,y1 (0)y2 (0)0,y1 (0)y2 (0)1,t0,0.2.運(yùn)用 3 中代碼，我們?nèi)〔介L(zhǎng)h=0.02, 我們求解得tx(t)x(t)y(t)y(t)0.0000000.0000001.0000000.0000001.0000000.0200000.0185090.7847200.0185090.7847200.0400000.0290490.2327450.0290490.2327450.060000

17、0.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. 解洛倫茲方程方程如下，使用不同的初始值，觀察差別dx10x10 y

18、,dtdy28 xy xz,dtdz8xy.dt3 z分別使用初值 x0 , y0 , z0 5,5,5, x0 , y0 , z0 5.001,5,5解查看附件我們查看初始值和結(jié)束值。 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

19、-12.89172411.71253449.993.5494584.58850818.66144-10.450006-18.17170016.666380150.004.3004856.27903317.322649-14.303620-21.25238326.374359我們發(fā)現(xiàn)盡管初值僅有很小的區(qū)別，但是終值有的巨大的差別（與0.001 相比）。.初值為 x0 , y0 , z0 5,5,5畫出 yz,xz,xy 有，.初值為 x0 , y0 , z0 5.001,5,5畫出 yz,xz,xy 有，.V. 邊值問題我們考慮這樣一個(gè)方程.y (t1)y 2y(1t 2 )e t , 0t1y

20、(0)1; y(1)3.14使用 3 中代碼求解有詳細(xì)答案參看附件。我們看看幾個(gè)均分點(diǎn)的值ty(t) 打靶法y(t) 差分法0.01.0000001.0000000.11.2392241.2392240.31.7030171.7030170.52.1442612.1442610.72.5609032.5609040.92.9528452.9528451.03.1400003.140000在算法上，打靶法計(jì)算量明顯高于差分法，但是打靶法具有更高的普適性。在進(jìn)行，有解析解的方程求解時(shí)，發(fā)現(xiàn)在步長(zhǎng)相同時(shí)，差分法具有更高的精度。畫出解的圖有，Shooting 解法.有限差分解法End.代碼：文件 ma

21、in_ode.cpp#include #include #include #include ode.h/#include ./array.hvoid free2DArray(double *a, int m)int i;for( i=0; im; i+ )free(ai);free(a);/ y= f(t,y), fun = f(t,y)double fun(double t,double y)return exp(t)/y;.void funv1(double t,double y,double fv)/ fv0= y0 + 2.*y1; /fv1=3*y0 + 2.*y1;/ Lotka-

22、Volterra equationfv0= y0 - y0*y1 - 0.1 *y0*y0;fv1= y0*y1 - y1 - 0.05*y1*y1;void funv2(double t,double y,double fv)/ y=x*x+x+1 fv0= y1;fv1= 4.*y0 - 4.*t*t - 4.*t - 2.;void funv3(double t, double y, double fv)fv0 = y1;fv1 = -278.28 / 0.15*y0 + 110.57 / 0.15*y2;fv2 = y3;fv3 = -278.28 / 0.15*y2 + 110.57

23、 / 0.15*y0;void funv4(double t, double y, double fv)fv0 = y1;fv1 = -(t + 1.)*y1 + y0 + (1. - t*t)*exp(-t);double pfun(double t)return -(t+1.);double qfun(double t)return 1.;double rfun(double t)return (1. - t*t)*exp(-t);/ -4.*t*t - 4.*t - 2.;void funvLorenz(double t,double yv,double fv)double x=yv0,

24、 y=yv1, z=yv2;fv0= 10.*(-x+y);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 t0=0.,y0=2., tn=5., *yv,*yv2, *yMa, *yFD, yv02=2.,1., yvLorenz3=5.,5.,5.; double yv34 = 0., 1., 0., 1. ;/ exact solution: y2 = 2*exp(x)+2/ 1st order OD

25、EM=21;yv = RungeKutta_Heum( fun, t0, y0, tn, M);printf(y(5.)=%16.12f, %16.12f n, yvM-1, fabs(yvM-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, yv2M-1, fabs(yv2M-1-sqrt(2.*exp(5.)+2.);free(yv2);/ ODE systemyMa = RKSystem4th( funv1, 2, 0.

26、, yv0, 30., 1000);/*yv00 = 21.;yv01 = -9.;yMa = RKSystem4th( funv2, 2, -5., yv0, 5., 3001);*/printf( y130=%f, y230=%f n, yMa9990,yMa9991);/*printf(erry1=%f,erry2=%f, 4. * exp(4.) + 2. * exp(-1.) - yMa990, 6. * exp(4.) - 2.*exp(-1.) -yMa991);printf(erry1=%f,erry2=%f, 31 - yMa30000, 11 - yMa30001);*/

27、free2DArray(yMa,100);yMa = RKSystem4th( funv1, 2, 0., yv0, 30., 1000);fp=fopen(lv.dat,w+);for(i=0;i1000;i+)fprintf(fp,%ft %fn,yMai0,yMai1);fclose(fp);free2DArray(yMa,1000);yMa = RKSystem4th(funv3, 4, 0., yv3, 0.2, 11);fp = fopen(fv3.dat, w+);for (i = 0; i11; i+)fprintf(fp,%ft%ft%ft%ft%fn,0.02*i,yMai

28、0,yMai1,yMai2,yMai3);fclose(fp);.free2DArray(yMa, 11);/ Lorenz equM= 1001;yMa = RKSystem4th( funvLorenz, 3, 0., yvLorenz, 50., M);fp = fopen(Lorenz01.dat,w+);for(i=0;iM;i+)fprintf(fp,%ft %ft %fn, yMai0,yMai1,yMai2);fclose(fp);M = 1001;yvLorenz0 = 5.001;yMa = RKSystem4th(funvLorenz, 3, 0., yvLorenz,

29、50., M);fp = fopen(Lorenz02.dat, w+);for (i = 0; iM; i+)fprintf(fp, %ft %ft %fn, yMai0, yMai1, yMai2); fclose(fp);/ high order ODE/ BVPM = 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;iM;i+)fprintf(fp, %ft %fn, t0 + (tn - t0) / (M-1) *

30、i, yMai0);fclose(fp);free2DArray(yMa,M);/BVP FDM = 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; iM; i+)fprintf(fp, %ft %fn, t0 +(tn-t0)/(M-1)*i, yFDi);fclose(fp);free(yFD);return 0;.文件 ode.cpp#include #include ode.h#include #

31、include /#include ./array.hdouble* RungeKutta_Heum( double fun(double,double), double t0, double y0, double tn, int M)/y(t+h)=y(t)+h/ 2*(f(t,y)+f(t+h,y+hf(t,y)double h = (tn - t0) / (M-1);double t = t0;/double y100 = 0 ;double *y;y = (double *)malloc(M * sizeof(double);int i = 0;y0 = y0;for (i = 0;

32、i M-1; i+)yi+1 = yi + h / 2 * (fun(t, yi) + fun(t + h, yi + h*fun(t, yi);t = t + h;return y;/* double* RungeKutta_EulerCauchy( double fun(double,double), double a, double b, double y0, int M) */double* RungeKutta4th( double fun(double,double), double t0, double y0, double tn, int M)double h = (tn -

33、t0) / (M - 1);double t = t0;/double y100 = 0 ;double *y;y = (double *)malloc(M * sizeof(double);int i = 0;y0 = y0;for (i = 0; i M-1; i+)yi + 1 = yi + h / 6 * (fun(t, yi) + 2 * fun(t + h / 2, yi + h / 2 * fun(t, yi) + 2 * fun(t + h / 2, yi + h / 2 * fun(t + h / 2, yi + h / 2 * fun(t, yi) + fun(t + h,

34、 yi+h*fun(t + h / 2, yi + h / 2 * fun(t + h / 2, yi + h / 2 * fun(t,yi);t = t + h;.return y;double* RKSystem4th( void fun(double,double,double ), int N, double t0, double y00,double tn, int M)double h = (tn - t0) / (M - 1);double t = t0;/double y100 = 0 ;double *y;double *dy;double *tempdy;y = (double *