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

1、使用C语言解常微分方程 C解常微分方程姓名：Vincent年级：2010，学号：1033*，组号：5（小组），4（大组）1. 数值方法:我们的实验目标是解常微分方程，其中包括几类问题。一阶常微分初值问题，高阶常微分初值问题，常微分方程组初值问题，二阶常微分方程边值问题，二阶线性常微分方程边值问题。对待上面的几类问题，我们分别使用不同的方法。 初值问题使用 龙格-库塔 来处理 边值问题用打靶法来处理 线性边值问题有限差分法初值问题我们分别使用 二阶 龙格-库塔 方法 4阶 龙格-库塔 方法来处理一阶常微分方程。理论如下：对于这样一个方程当h很小时，我们用泰勒展开，当我们选择正确的参数 aij,b

2、ij之后，就可以近似认为这就是泰勒展开。对于二阶，我们有：其中经过前人的计算机经验，我们有，选择 A=1/2,B=1/2,则 P=1,Q=1，于是又如下形式，这也使休恩方法的表达式。所以我们称其为 龙格（库塔）休恩方法对于4阶龙格方法，我们有类似的想法，我们使用前人经验的出的系数，有如下公式对于高阶微分方程及微分方程组我们用 4阶龙格-库塔方法来解对于一个如下的微分方程组我们可以认为是一个一阶向量微分方程，所以可以用龙格-库塔方法的向量形式解。对于一个高阶的微分方程，形式如下:我们可以构建出一个一阶的微分方程组，令则有所以我们实际只要解一个微分方程组 其中初值为 使用4阶龙格-库塔方法， 使用

3、这个向量形式的龙格-库塔方法我们便可就出方程的数值解。边值问题对于边值问题，我们分为两类 一般的边值问题 线性边值微分方程一般的边值问题，我们是使用打靶法来求解，对于这样一个方程主要思路是，化成初值问题来求解。我们已有 这样我们便可求出方程的解。 线性微分方程边值问题对于这样的问题，我们可以使用一个更加高效的方法，有限差分法。对于如下方程 我们对其进行差分这样的话，我们的微分方程可以写成，于是我们得到了个线性方程组这样的话我们求解出x对于上面的矩阵，我们可以分解出两个三角阵，然后回代求解便可。我们用回代法可以直接求解至此，我们便求出了目标方程的解2. 流程图 二阶龙格-库塔对于i = 0到M-

4、1; yi+1 = yi + h / 2 * (fun(t, yi) + fun(t + h, yi + h*fun(t, yi);return y; 4阶龙格-库塔对于i = 0到M-1;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, yi+h*fun(t + h / 2, yi + h / 2 * fun(t +

5、h / 2, yi + h / 2 *fun(t, yi);return y; 4阶龙格-库塔解方程组对于k= 0到M-1; 对于i= 0到N; fun(t, yk, dy0) 对于i= 0到N; tempdyj = ykj + h / 2 * dy0j; fun(t + h / 2, tempdy, dy1); 对于i= 0到N; tempdyj = ykj + h / 2 * dy1j; fun(t + h / 2, tempdy, dy2); 对于i= 0到N; tempdyj = ykj + h * dy2j; fun(t + h, tempdy, dy3); yk+1i = yki

6、+ h / 6 * (dy0i + 2 * dy1i + 2 * dy2i + dy3i); return y; 打靶法当errepsilon y = 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 = fabs(u1 - v1); return y; 有限差分法对i 从 0到 m；求

7、出，b,r,a,c bi = 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) - 1.)*alpha; rm - 1 = rm - 1 - (h / 2 * pfun(tn - h) - 1)*beta;求出d,ud0 = b0; u0 = c0 / b0;对i 从 1到m - 1 di = bi - ai * ui - 1; ui = ci / di;dm - 1 = bm - 1 - am - 1

8、 * 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. 解下面的方程求t=5时，y的值 使用3中代码执行得到My(5) 2阶方法解y(5)4阶方法解2阶方法误差4阶方法误差100.204080 0. 00050. 000000214374160对比发现4阶精度远高于2阶精度当我们细分到一定程度之后，我们发

9、现误差的减小慢慢变小。所以，若需要极高的精度的话会花费大量的计算资源。II. 解下面的方程组我们选择M=1000来细分运用3中代码，我们求解得 III. 解下面高阶微分方程(震动方程)运用3中代码，我们取步长h=, 我们求解得画出解y1,y2有，IV. 解洛伦兹方程方程如下，使用不同的初始值，观察差别分别使用初值解查看附件我们查看初始值和结束值。txyzxyz055555我们发现尽管初值仅有很小的区别，但是终值有的巨大的差别（与相比）。初值为画出yz,xz,xy有，初值为画出yz,xz,xy有，V. 边值问题我们考虑这样一个方程使用3中代码求解有详细答案参看附件。我们看看几个均分点的值ty(t

10、) 打靶法y(t)差分法在算法上，打靶法计算量明显高于差分法，但是打靶法具有更高的普适性。在进行，有解析解的方程求解时，发现在步长相同时，差分法具有更高的精度。画出解的图有，Shooting 解法 有限差分解法End.代码：文件 #include #include #include #include /void free2DArray(double *a, int m) int i; for( i=0; im; i+ ) free(ai); free(a);y1; y1; y0 - 4.*t*t - 4.*t - 2.;void funv3(double t, double y, double

11、 fv) fv0 = y1; fv1 = / *y0 + / *y2; fv2 = y3; fv3 = / *y2 + / *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);t*t - 4.*t - 2.;voi

12、d funvLorenz(double t,double yv,double fv) double x=yv0, 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.

13、, 1., 0., 1. ; =%, % 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.)=%, % n, yv2M-1, fabs(yv2M-1-sqrt(2.*exp(5.)+2.); free(yv2); yv0, 30., 1000);/* yv00 = 21.; yv01 = -9.; yMa = RKSystem4th( funv2, 2, -5., yv0, 5., 3001);*/ printf( y130

14、=%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);*/ free2DArray(yMa,100); yMa = RKSystem4th( funv1, 2, 0., yv0, 30., 1000); fp=fopen(,w+); for(i=0;i1000;i

15、+) fprintf(fp,%ft %fn,yMai0,yMai1); fclose(fp); free2DArray(yMa,1000); yMa = RKSystem4th(funv3, 4, 0., yv3, , 11); fp = fopen(, w+); for (i = 0; i11; i+) fprintf(fp, %ft %ft %ft %ft %fn,*i, yMai0, yMai1, yMai2, yMai3); fclose(fp); free2DArray(yMa, 11); yvLorenz, 50., M); fp = fopen(,w+); for(i=0;iM;

16、i+) fprintf(fp,%ft %ft %fn, yMai0,yMai1,yMai2); fclose(fp); M = 1001; yvLorenz0 = ; yMa = RKSystem4th(funvLorenz, 3, 0., yvLorenz, 50., M); fp = fopen(, w+); for (i = 0; iM; i+) fprintf(fp, %ft %ft %fn, yMai0, yMai1, yMai2); fclose(fp); tn = 1.; a = 1.; b = ; yMa = BVP_Shooting(funv4, t0, a, tn, b,

17、M); fp=fopen(,w+); for(i=0;iM;i+) fprintf(fp, %ft %fn, t0 + (tn - t0) / (M-1) *i, yMai0); 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; iM; i+) fprintf(fp, %ft %fn, t0 +(tn-t0)/(M-1)*i, yFDi); fclose(fp); free(yF

18、D); return 0;文件#include #include #include #include /double* RungeKutta_Heum( double fun(double,double), double t0, double y0, double tn, int M) ; double v2 = a, -2.; double *y; double x2,err; double f0,f1; do y = RKSystem4th( fun, 2, t0, u, tn, M); f0 = yM-10 - b; u1 = y01; y = RKSystem4th( fun, 2,

19、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); i+; while( err epsilon); return y;/*double* BVP_FD( double pfun(double), double qfun(double), double rfun(double), double t00,double alpha, double tn,double beta, int m ) - h * h *qfun(t);

20、ri =h * h * rfun(t); t = t + h; r0 = r0 - (h/2. * pfun(t0) + 1.) * alpha; - h/2. * pfun(tn) * beta; t = t0; for (i = 0; i M - 1; i+) ai = 1. + h / 2. * pfun(t + h); ci = 1. - h / 2. * pfun(t); *( pfun(t0+h) - 1.)*alpha; rm - 1 = rm - 1 - (h / 2 * pfun(tn - h) - 1)*beta; d0 = b0; u0 = c0 / b0; for (i

21、 = 1; i m - 1; i+) di = bi - ai * ui - 1; if (di = 0) return 0; ui = ci / di; dm - 1 = bm - 1 - am - 1 * um - 2; y0 = r0 / d0; for (i = 1; i = 0;i- ) xi + 1 = yi - ui * xi + 2; x0 = alpha; xM - 1 = beta; return x;double* BVP_Lshooting( double pfun(double), double qfun(double), double rfun(double), d

22、ouble t0,double a, double tn,double b, int M ) return 0;文件#ifndef ODE_HHHH#define ODE_HHHH/ 2nd order Runge-Kutta method for solving initial value problemdouble* RungeKutta_Heum( double fun(double,double), double t0, double y0, double tn, int M);/ 4th order Runge-Kutta method for solving initial val

23、ue problemdouble* RungeKutta4th( double fun(double,double), double t0, double y0, double tn, int M);/ ODE systems and high order ODE solverdouble* RKSystem4th( void fun(double,double,double ), int N, double t0, double y0,double tn, int M);/ ODE boundary value problem using shooting methoddouble* BVP_Shooting( void fun(double,double,double), double t0,double a, double tn,double b, int M );/ BVP using finite difference methoddouble* BVP_FD( double pfun(double), double qfun(double), double rfun(double), double t0,double a, double tn,double b, int M );#endif