使用C语言解常微分方程 C ODE

《使用C语言解常微分方程 C ODE》由会员分享，可在线阅读，更多相关《使用C语言解常微分方程 C ODE（15页珍藏版）》请在装配图网上搜索。

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

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

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

4、) + 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 + h / 2, yi + h / 2 *fun(t, yi);ret

5、urn 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 + h / 6 * (dy0i + 2 * dy1i + 2 * dy2i + dy3i)

6、;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;*2= u1 - f0 * (v1-u1)/(f1-f0);u1 = v1;v1 = *2;err = fabs(u1 - v1);return y; 有限差分法对i 从 0到 m；求出，b,r,a,cbi = 2 + h*h*qfun(t);ri = -h*h*rfun(t);ai = -h /

7、 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 - 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;*m

8、 = ym -1;对i 从 m 2到0*i + 1 = yi - ui * *i + 2;*0 = alpha;*M - 1 = beta;return *;3. 代码实现过程查看4. 数值例子与分析I. 解下面的方程求t=5时，y的值使用3中代码执行得到My(5) 2阶方法解y(5)4阶方法解2阶方法误差4阶方法误差1017.717.10.70.02017.517.70.40.74017.017.30.90.28017.6 17.5 0. 50. 416017.917.20.80.2比照发现4阶精度远高于2阶精度当我们细分到一定程度之后，我们发现误差的减小慢慢变小。所以，假设需要极高的精度的

9、话会花费大量的计算资源。II. 解下面的方程组我们选择M=1000来细分运用3中代码，我们求解得III. 解下面高阶微分方程(震动方程)运用3中代码，我们取步长h=0.02, 我们求解得画出解y1,y2有，IV. 解洛伦兹方程方程如下，使用不同的初始值，观察差异分别使用初值解查看我们查看初始值和完毕值。t*yz*yz05555.001550.016.2021 960.6576426.3873596.20241310.6585266.3877950.029.37401417.4526590.7163959.37497817.45400710.71778349.983.2582763.369219

10、20.608133-7.008305-12.89172411.71253449.993.5494584.58850818.661441-10.450006-18.17170016.66638050.004.3004856.27903317.322649-14.303620-21.25238326.374359我们发现尽管初值仅有很小的区别，但是终值有的巨大的差异与0.001相比。初值为画出yz,*z,*y有，初值为画出yz,*z,*y有，V. 边值问题我们考虑这样一个方程使用3中代码求解有详细答案参看。我们看看几个均分点的值ty(t) 打靶法y(t)差分法0.01.0000001.000000

11、0.11.2392241.2392240.31.7030171.7030170.52.1442612.1442610.72.5609032.5609040.92.9528452.9528451.03.1400003.140000在算法上，打靶法计算量明显高于差分法，但是打靶法具有更高的普适性。在进展，有解析解的方程求解时，发现在步长一样时，差分法具有更高的精度。画出解的图有，Shooting 解法有限差分解法End.代码：文件 main_ode.cpp*include *include *include *include ode.h/*include ./array.hvoid free2DA

12、rray(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 e*p(t)/y;void funv1(double t,double y,double fv)/ / fv0= y0 + 2.*y1; /fv1=3*y0 + 2.*y1;/ Lotka-Volterra equationfv0= y0 - y0*y1 - 0.1 *y0*y0;fv1= y0*y1 - y1 - 0.05*y1*y1;void f

13、unv2(double t,double y,double fv)/ y=*+*+1fv0= 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 / 0.15*y0;void funv4(double t, double y, double fv)fv0 = y1;fv1 = -(t + 1.)*y1 + y0

14、 + (1. - t*t)*e*p(-t);double pfun(double t)return -(t+1.);double qfun(double t)return 1.;double rfun(double t)return (1. - t*t)*e*p(-t);/ -4.*t*t - 4.*t - 2.;void funvLorenz(double t,double yv,double fv)double *=yv0, y=yv1, z=yv2;fv0= 10.*(-*+y);fv1= 28.* - y - *z;fv2= -8./3.*z + *y;int main( int ar

15、gc, 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. ;/ e*act solution: y2 = 2*e*p(*)+2/ 1st order ODEM=21;yv = RungeKutta_Heum( fun, t0, y0, tn, M);printf(y(5.)=%16.12f, %16.12f n, yvM-1, f

16、abs(yvM-1-sqrt(2.*e*p(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.*e*p(5.)+2.);free(yv2);/ ODE systemyMa = RKSystem4th( funv1, 2, 0., yv0, 30., 1000);/*yv00 = 21.;yv01 = -9.;yMa = RKSystem4th( funv2, 2, -5., yv0, 5., 3001)

17、;*/printf( y130=%f, y230=%f n, yMa9990,yMa9991);/*printf(erry1=%f,erry2=%f, 4. * e*p(4.) + 2. * e*p(-1.) - yMa990, 6. * e*p(4.) - 2.*e*p(-1.) - yMa991);printf(erry1=%f,erry2=%f, 31 - yMa30000, 11 - yMa30001);*/free2DArray(yMa,100);yMa = RKSystem4th( funv1, 2, 0., yv0, 30., 1000);fp=fopen(lv.dat,w+);

18、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, yMai0, yMai1, yMai2, yMai3);fclose(fp);free2DArray(yMa, 11);/ Lorenz equM = 1001;yMa = R

19、KSystem4th( 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, 50., M);fp = fopen(Lorenz02.dat, w+);for (i = 0; iM; i+)fprintf(fp, %ft %ft %fn,

20、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) *i, yMai0);fclose(fp);free2DArray(yMa,M);/BVP FDM = 101;t0 = 0.;tn = 1.;a = 1.;b =

21、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 *include /*include ./array.hdouble* RungeKutta_Heum( double fun(double,double), doubl

22、e 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; i M-1; i+)yi+1 = yi + h / 2 * (fun(t, yi) + fun(t + h, yi + h*fun(t, yi);t = t + h;re

23、turn 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 - t0) / (M - 1);double t = t0;/double y100 = 0 ;double *y;y = (double *)malloc(M * size

24、of(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, yi+h*fun(t + h / 2, yi + h / 2 * fun(t + h / 2, yi + h / 2 * fun(t, yi);t = t + h;ret

25、urn 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 *)malloc(M * sizeof(double*);dy = (double *)malloc(4 * sizeof(double*);tempdy = (do

26、uble *)malloc(N*sizeof(double);int i = 0, j = 0, k = 0;for (i = 0; i M; i+)yi = (double *)malloc(N * sizeof(double);for (i = 0; i 4; i+)dyi = (double *)malloc(N*sizeof(double);for (i = 0; i N; i+)y0i = y00i;for (k = 0; k M-1; k+)for (i = 0; i N; i+)fun(t, yk, dy0);for (j = 0; j N; j+)tempdyj = ykj +

27、 h / 2 * dy0j;fun(t + h / 2, tempdy, dy1);for (j = 0; j N; j+)tempdyj = ykj + h / 2 * dy1j;fun(t + h / 2, tempdy, dy2);for (j = 0; j 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 )/y=p(t)y+q(t)

28、y+r(t)int i = 0;int M = m - 2;double h = (tn - t00) / (M+1);double t0 = t00 + h;double t = t0;/(-h/2*pj-1)*j-1+(2+h*h*qj)*j+(h/2*pj-1)*j+1=-h*h*rj;double *a, *b, *c, *d, *u, *r, *y, *; /double *l;a = (double *)malloc(M - 1) * sizeof(double);b = (double *)malloc(M * sizeof(double);c = (double *)mallo

29、c(M - 1) * sizeof(double);/l = (double *)malloc(M - 1) * sizeof(double);d = (double *)malloc(M * sizeof(double);u = (double *)malloc(M - 1) * sizeof(double);r = (double *)malloc(M * sizeof(double);y = (double *)malloc(M * sizeof(double);* = (double *)malloc(M * sizeof(double);t = t0;for(i = 0; i M;

30、i+)bi = -2. - h * h *qfun(t);ri =h * h * rfun(t);t = t + h;r0 = r0 - (h/2. * pfun(t0) + 1.) * alpha;/r0 = 1.0855;rM - 1 = rM - 1 - (1. - 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);/li = ai;t = t + h;d0 = b0;u0 = c0 / b0;for (i = 0;

31、 i M - 2; i+)di + 1 = bi + 1 - ai * ui;if (0 = di + 1)printf(error!n);return 0;ui + 1 = ci + 1 / di + 1;dM - 1 =bM - 1-aM - 2 * uM - 2;/回代y0 = r0 / d0;for (i = 1; i = 0; i-)*i = yi - ui * yi + 1;return *;*/double* BVP_FD(double pfun(double), double qfun(double), double rfun(double),double t0, double

32、 alpha, double tn, double beta, int M)int i = 0;double h = (tn - t0) / (M-1);double t = t0;int m = M - 2;double *a, *b, *c, *d, *u, *r, *y, *;a = (double *)malloc(m * sizeof(double);b = (double *)malloc(m * sizeof(double);c = (double *)malloc(m * sizeof(double);d = (double *)malloc(m * sizeof(double

33、);u = (double *)malloc(m * sizeof(double);r = (double *)malloc(m * sizeof(double);y = (double *)malloc(m * sizeof(double);* = (double *)malloc(M * sizeof(double);t = t0 + h;for (i = 0; i m; i+)bi = 2 + h*h*qfun(t);ri = -h*h*rfun(t);ai = -h / 2 * pfun(t) - 1;ci = h / 2 * pfun(t) - 1;t = t + h;r0 = r0

34、 - (-h / 2. )*( pfun(t0+h) - 1.)*alpha;rm - 1 = rm - 1 - (h / 2 * pfun(tn - h) - 1)*beta;d0 = b0;u0 = c0 / b0;for (i = 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- )*i + 1 = yi - ui * *i + 2;*0 = alpha;*M -

35、1 = beta;return *;double* BVP_Lshooting( double pfun(double), double qfun(double), double rfun(double), double t0,double a, double tn,double b, int M )return 0;文件ode.h*ifndef ODE_HHHH*define ODE_HHHH/ 2nd order Runge-Kutta method for solving initial value problemdouble* RungeKutta_Heum( double fun(d

36、ouble,double), double t0, double y0, double tn, int M);/ 4th order Runge-Kutta method for solving initial value 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 )

37、, 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. z.