使用C语言解常微分方程 C ODEWord文件下载.docx

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

2、程的解2. 流程图 二阶龙格-库塔对于i = 0到M-1; yi+1 = yi + h / 2 * （fun（t, yi） + fun（t + h, yi + h*fun（t, yi）;return y; 4阶龙格-库塔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

3、 / 2 * fun（t + h / 2, yi + h / 2 *fun（t, yi）; 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）; tempdyj = ykj + h * dy2j; fun（t + h, tempdy, dy3）; yk+1i = yki + h

4、/ 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）; 有限差分法对i 从 0到 m；求出，b,r,a,c bi =

5、 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 * um - 2;回代y0

6、 = 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阶方法误差1017.236717.64310.48670.00002017.997517.08970.04744017.145017.13230.00.000003443822

7、8017.5806 17.1875 0. 00050. 00000021437416017.163917.08720.00080.000000013372对比发现4阶精度远高于2阶精度当我们细分到一定程度之后，我们发现误差的减小慢慢变小。所以，若需要极高的精度的话会花费大量的计算资源。II. 解下面的方程组我们选择M=1000来细分运用3中代码，我们求解得III. 解下面高阶微分方程（震动方程）运用3中代码，我们取步长h=0.02, 我们求解得画出解y1,y2有，IV. 解洛伦兹方程方程如下，使用不同的初始值，观察差别分别使用初值解查看附件我们查看初始值和结束值。txyz55.0010.016

8、.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.661441-10.450006-18.17170016.66638050.004.3004856.27903317.322649-14.303620-21.25238326.374359我们发现尽管初值仅有很小的区别，但是终值有

9、的巨大的差别（与0.001相比）。初值为画出yz,xz,xy有，V. 边值问题我们考虑这样一个方程使用3中代码求解有详细答案参看附件。我们看看几个均分点的值y（t） 打靶法y（t）差分法1.0000000.11.2392240.31.7030170.52.1442610.72.5609032.5609040.92.9528451.03.140000在算法上，打靶法计算量明显高于差分法，但是打靶法具有更高的普适性。在进行，有解析解的方程求解时，发现在步长相同时，差分法具有更高的精度。画出解的图有，Shooting 解法有限差分解法End.代码：文件 main_ode.cpp#include st

10、dlib.hmath.h#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-Volterra equation

11、 fv0= 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 / 0.15*y

12、0;void funv4（double t, double y, double fv） 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, y=yv1, z=yv2

13、; 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 ODE M=21

14、; 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）; yv2= RungeKutta4th（ fun, t0, y0, tn, M）;, yv2M-1, fabs（yv2M-1-sqrt（2.*exp（5.）+2.）; free（yv2）; / ODE system yMa = RKSystem4th（ funv1, 2, 0., yv0, 30., 1000）;/* yv00 = 21.; yv

15、01 = -9.; yMa = RKSystem4th（ funv2, 2, -5., yv0, 5., 3001）;*/ 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）;, 31 - yMa30000, 11 - yMa30001）; free2DArray（yMa,100）; fp=fopen（lv.dat,w+）; for（i=0;i1000;i+） fp

16、rintf（fp,%ft %fn,yMai0,yMai1）; fclose（fp）; free2DArray（yMa,1000）; yMa = RKSystem4th（funv3, 4, 0., yv3, 0.2, 11）; fp = fopen（fv3.dat, for （i = 0;11; i+） fprintf（fp, %ft %ft %ft %ft %fn,0.02*i, yMai0, yMai1, yMai2, yMai3）; free2DArray（yMa, 11）; / Lorenz equ M = 1001; yMa = RKSystem4th（ funvLorenz, 3,

17、0., yvLorenz, 50., M）;Lorenz01.datM;%ft %ft %fn, yMai0,yMai1,yMai2）; yvLorenz0 = 5.001; yMa = RKSystem4th（funvLorenz, 3, 0., yvLorenz, 50., M）;Lorenz02.dat, yMai0, yMai1, yMai2）; / high order ODE / BVP M = 101; t0 = 0.; tn = 1.; a = 1.; b = 3.14; yMa = BVP_Shooting（funv4, t0, a, tn, b, M）;Shoot.dat,

18、 t0 + （tn - t0） / （M-1） *i, yMai0）; free2DArray（yMa,M）; /BVP FD yFD = BVP_FD（pfun, qfun, rfun, t0, a, tn, b, M）;yFD.dat, t0 +（tn-t0）/（M-1）*i, yFDi）; free（yFD）; return 0;文件ode.cppdouble* 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

19、,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; i M-1; yi+1 = yi + h / 2 * （fun（t, yi） + fun（t + h, yi + h*fun（t, yi）; t = t + h; /* double* RungeKutta_EulerCauchy（ double fun（double,double）, double a,

20、 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）; 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,

21、 yi） + fun（t + h, yi+h*fun（t + h / 2, yi + h / 2 * fun（t + h / 2, yi + h / 2 * fun（t, yi）;double* RKSystem4th（ void fun（double,double,double ）, int N, double t0, double y00,double tn, int M） double *y; double *dy; double *tempdy; y = （double *）malloc（M * sizeof（double*）; dy = （double *）malloc（4 * sizeof（double*）; tempdy = （double *）malloc（N*sizeof（double）; int i = 0, j = 0, k = 0; M; yi = （double *）malloc（N * sizeof（double）; 4; dyi = （double *）malloc（N*sizeof（double）; N; y0i = y00i; for （k = 0; k k+） for （i = 0; fun（t, yk, dy0）; for （j = 0; j j+） tempdyj = y