第8章常微分方程数值解法

第8章常微分方程数值解法8.1引言在常微分方程中利用典型方程的解法，对这些方程就可以求得其解的解析表达式。例如，对一阶线性微分方程(8.1)有通解为有了解析解，就可以根据初始条件或边界条件把其中的任意常数具体确定下来。)()(xvyxuycx)x(vyxd)x(uxd)x(udee在许多工程实践和科学研究中，除了少数特殊类型的微分方程能用解析法求得其精确解外，大多数情况要得出解的解析表达式是极其困难的，甚至是不可能的，因此，只能用近似方法求得其近似解。例如我们很容易求出初值问题的解为但要计算它的值，还需要用数值积分的方法。如果要对许多个值计算解的近似值，那么工作量非常大。况且实际计算不一定要求解析表达式，而是只需求在某些点上满足精度的解的近似值或解的近似表达式就可以了。0(0)21yxyyxtxtxy0dee)(22x)(xy由于高阶常微分方程可以转化为一阶常微分方程组，因此，为了不失一般性，本章主要介绍一类一阶常微分方程初值问题在区间上的数值计算方法。在初值问题的数值解中，常见的有单步法和多步法。所谓单步法，就是在算法中根据初始条件，由解就能求出；所谓多步法，就是在求时不仅要知道，而且还要知道在前面的值。本章首先介绍单步法的欧拉法及龙格-库塔法，然后介绍多步法的阿当姆斯法，它们是解初值问题的几个最基本的常用的有效方法。bxa0yiy)1,2,(1iyi1iyiyiy,y,y2i1i00)()(yxyyx,fy8.2欧拉（Euler）法8.2.1欧拉法的基本思想虽然欧拉法在实际应用时较少采用，但由于它的结构简单，在理论上具有非常重要的意义。欧拉法的基本思想就是用差分方程初值问题(8.3)的解来近似微分方程初值问题(8.2)的解，其中，式(8.3)也称为欧拉公式。)()0,1,()(01ayyiy,xfhyyiiii2)-(/abh欧拉法的几何意义是用一条自点出发的折线去逼近积分曲线，如图8.1所示。因此，这种方法又称为折线法。显然，欧拉法简单地取折线的端点作为数值解，精度非常差。)(00y,x)(xfy8.2.2实现欧拉法的基本步骤（1）输入的区间，区间的等分数以及在处的初始值；（2）对，计算：（3）输出；（4）结束。xnx,x0Ny0x0y10,1,2,Ni)(iiiy,xffiiifhyyhi*xxiiN/xxhn)-(0Ny例8.1用欧拉法求解初值问题的数值解（保留4位小数）。02(0)y10.y,xxyy#includestdio.hfloatfunc(floatx,floaty){return(y-x);}floateuler(floatx0,floatxn,floaty0,intN){floatx,y,h;inti;x=x0;y=y0;h=(xn-x0)/(float)N;/*计算步长*/for(i=1;i=N;i++)/*欧拉公式*/{y=y+h*func(x,y);x=x0+i*h;}return(y);}main(){floatx0,xn,y0,e;intN;clrscr();printf(\ninputn:\n);/*输入区间等分数*/scanf(%d,&N);printf(inputx0,xn:\n);/*输入x的区间*/scanf(%f%f,&x0,&xn);printf(inputy0:\n);/*输入x0处的y的值*/scanf(%f,&y0);e=euler(x0,xn,y0,N);printf(y(%f)=%6.4f,y0,e);}程序运行结果：inputn:5inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.4883inputn:10inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.5937inputn:20inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.65331e)(xxyx此问题的精确解为，相应计算值为4.7183。从计算结果可以看出，虽然等分数越大，数值解越精确，但精度并不理想。8.3改进的欧拉法8.3.1改进的欧拉法的基本思想改进的欧拉法以隐函数的方式给出，计算时常使用迭代法，一般可先由欧拉公式(8.3)计算出的初始值，再按下面的格式进行的迭代计算(8.4)改进的欧拉法由于添加了迭代过程，计算量较欧拉法增加了一倍，付出这种代价的目的是为了提高精度。1iy1iy(0)1iy1iy)0,1,2,()()(2)()(111)(1(0)1ky,xfy,xfhyyy,xfhyykiiiiikiiiii8.3.2实现改进的欧拉法的基本步骤（1）输入的区间，区间的等分数以及在处的初始值；（2）,，；（3）；（4）计算：（5）若，，转（4）；否则，转（6）；（6）输出；（7）结束。xnx,x0Ny0x0y/N)x-x(h0n0xx0yy1i)yx,(fhyyphi*xx0)(pcyx,fh*yy2(/)yyycpNi1iiy例8.2用改进的欧拉法求解初值问题的数值解（保留4位小数）。02(0)y10.y,xxyy#includestdio.hfloatfunc(floatx,floaty){return(y-x);}floateuler(floatx0,floatxn,floaty0,intN){floatx,y,yp,yc,h;inti;x=x0;y=y0;h=(xn-x0)/(float)N;/*计算步长*/for(i=1;i=N;i++)/*改进的欧拉公式*/{yp=y+h*func(x,y);x=x0+i*h;yc=y+h*func(x,yp);y=(yp+yc)/2.0;}return(y);}main(){floatx0,xn,y0,e;intN;clrscr();printf(\ninputn:\n);/*输入区间等分数*/scanf(%d,&N);printf(inputx0,xn:\n);/*输入x的区间*/scanf(%f%f,&x0,&xn);printf(inputy0:\n);/*输入x0处的y的值*/scanf(%f,&y0);e=euler(x0,xn,y0,N);printf(y(%f)=%6.4f,y0,e);}程序运行结果：inputn:5inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.7027inputn:10inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.7141inputn:20inputx0,xn:0.01.0inputy0:2.0y(2.000000)=4.71728.4龙格—库塔（Runge-Kutta）法8.4.1龙格—库塔法的基本思想在欧拉法中，用解函数在点处的斜率计算从到的增量，的表达式与的Taylor展开式的前二项相等，使方法只有一阶精度。改进的欧拉法用两个点，处的斜率、的平均值计算增量，使方法具有二阶精度，即的表达式与的Taylor展开式的前三项相等。由此龙格和库塔提出了一种间接地运用Taylor公式的方法，即利用在若干个待定点上的函数值和导数值做出线性组合式，选取适当系数使这个组合式进行Taylor展开后与的Taylor展开式有较多的项达到一致，从而得出较高阶的数值公式，这就是龙格—库塔法的基本思想。)0,1,()(1iy,xfhyyiiii)(xfyix)(iiy,xfiy1iy1iy)(1ixyix1ix)(iiy,xf)(11iiy,xf1iy)(1ixy)(xy)(1ixy8.4.2实现标准四阶龙格—库塔法的基本步骤（1）输入的区间，区间的等分数以及在处的初始值；（2）,，；（3）；（4）计算：（5）若，，转（4）；否则，转（6）；（6）输出；（7）结束。xnx,x0Ny0x0y/Nxxhn)-(00xx0yy1i2/hxxh)(1yx,fd)21(2/dh*y,xfdh)22(3/dh*y,xfdh)3(4dh*y,xfdh4)/63*22*21(ddddh*yyhi*xx0Ni1iiy例8.3用标准四阶龙格—库塔法求解初值问题的数值解（保留6位小数）。0(0)10yxyxy#includestdio.hfloatfunc(floatx,floaty){return(x-y);}floatrunge_kutta(floatx0,floatxn,floaty0,intN){floatx,y,y1,y2,h,xh;floatd1,d2,d3,d4;inti;x=x0;y=y0;h=(xn-x0)/(float)N;/*计算步长*/for(i=1;i=N;i++)/*标准四阶龙格-库塔公式*/{xh=x+h/2;d1=func(x,y);d2=func(xh,y+h*d1/2.0);d3=func(xh,y+h*d2/2.0);d4=func(xh,y+h*d3);y=y+h*(d1+2*d2+2*d3+d4)/6.0;x=x0+i*h;}return(y);}main(){floatx0,xn,y0,e;intN;clrscr();printf(\ninputn:\n);/*输入区间等分数*/scanf(%d,&N);printf(inputx0,xn:\n);/*输入x的区间*/scanf(%f%f,&x0,&xn);printf(inputy0:\n);/*输入x0处的y的值*/scanf(%f,&y0);e=runge_kutta(x0,xn,y0,N);printf(y(%f)=%8.6f,y0,e);}程序运行结果：inputn:5inputx0,xn:0.01.0inputy0:0.0y(0.000000)=0.356261inputn:10inputx0,xn:0.01.0inputy0:0.0y(0.000000)=0.362344inputn:20inputx0,xn:0.01.0inputy0:0.0y(0.000000)=0.365179注意：此问题的精确解为0.367879。从计算结果可以看出，随着等分数的增加，数值解越精确。但需要指出的是：龙格—库塔法的推导是在用Taylor展开法的基础上进行的，因而它要求所求的解具有较好的光滑性质。假若解的光滑性差，那么使用标准四阶龙格—库塔法求得的数值解，其精度可能反而不如改进的欧拉法。因此，在实际计算时，应当针对问题的具体特点，选择合适的算法。N8.5阿当姆斯（Adams）法8.5.1阿当姆斯法的基本思想龙格—库塔法是一类重要算法，但是，这类算法在每一步都要先预报几个点上的斜率值，因而计算量比较大。事实上，考虑到在计算之前已得出一系列节点上的斜率值，能否利用这些已知信息来减少计算量，这正是阿当姆斯法的基本思想。构成阿当姆斯法的公式有阿当姆斯外插公式和阿当姆斯预测—校正公式两种，本节主要介绍阿当姆斯预测—校正公式。1iy,x,x,xiii218.5.2实现阿当姆斯法的基本步骤（1）输入的区间，区间的等分数以及在处的初始值；（2）；（3）使用标准四阶龙格—库塔法计算,为阿当姆斯预测—校正公式提供初始值；（4）（5）；（6）若，，转（5）；否则，转（7）；（7）输出；（8）结束。xnx,x0Ny0x0y/Nxxhn)-(0]2yy[]1yy[]0yy[,,]2yy[]1yy[2]0yy[1y,y,y3ihi*xx0))/2403(9-)12(37)2(59-)(55(yh,xfyh,xfyh,xfyx,fh*ypy))/2412()2(5-)(19)((9yh,xfyh,xfyx,fpyh,xfh*ycyc22110yy,yy,yy