补充3：常微分方程数值解法课件

数值计算方法第七章常微分方程数值解法1数值计算方法第七章常微分方程数值解法1在工程和科学计算中，所建立的各种常微分方程的初值或边值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因此，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究常微分方程的数值解法是有重要意义的。2在工程和科学计算中，所建立的各种常微分方程的7.0基本概念1.一阶常微分方程的初值问题(7.0-1)注：若f在D={a

x

b,|y|<+}内连续，且满足Lip条件：L

0，使

|f(x，y1)–f(x，y2)|

L|y1–y2|(7.0-2)则(7.0-1)的连续可微解y(x)在[a，b]上唯一存在。37.0基本概念32.初值问题的数值解称(7.0-1)的解y(x)在节点xi处的近似值yi

y(xi)a<x1<x2<...<xn=b.为其数值解，方法称为数值方法。注：①考虑等距节点:xi=a+ih，h=(b–a)/n.②从初始条件y(a)=y0出发，依次逐个计算y1，y2，…，yn的值，称为步进法。两种：单步法、多步法。42.初值问题的数值解4③二阶常微分方程y''(x)=f(x，y(x)，y'

(x))可设为一阶常微分方程组的初值问题：引进新的未知函数z(x)=y'(x)，则其初始条件为：称为一阶微分方程组的初值问题，方法类似。④边界问题，常用差分方法解。5③二阶常微分方程y''(x)=f(x，y(x)，y7.1初值问题数值解法的构造及其精度7.1.1构造方法对于(7.0-1)可借助Taylor展开（导数法）、差商法、积分法实现离散化来构造差分公式：1.Taylor展开：设y

C[a，b]将y(xi+1)=y(xi+h)在xi处展开

[xi，xi+1]

y(xi+1)

yi+hf(xi，yi)其中yi

y(xi).称yi+1=yi+hf(xi，yi).i=0，1，2，...，n–1(7.1-1)为Euler求解公式，（Euler法）67.1初值问题数值解法的构造及其精度62.用差商来表示：得差分方程：

yi+1=yi

+hf(xi，yi).即为Euler公式。若记

yi+1=yi

+hf(xi+1，yi+1).(7.1-2)称为向后Euler法。注：①Euler法为显式，向后Euler法为隐式——须解出yi+1.②可用迭代法yi+1(k+1)=yi

+hf(xi+1，yi+1(k))

k=0，1，2，…解得yi+1，其中yi+1(0)=yi+hf(xi，yi).72.用差商来表示：73.积分法:对(7.0-1)两边取积分得(7.1-3)取不同的数值积分可得不同的求解公式，如：①用矩形公式：左矩形

y(xi+1)

y(xi)+hf(xi，y(xi))Euler公式右矩形

y(xi+1)

y(xi)+hf(xi+1，y(xi+1))向后Euler公式83.积分法:对(7.0-1)两边取积分得8②用梯形公式：

(7.1-4)称(7.1-4)为梯形公式隐式公式。显化：预估值：校正值：(7.1-5)称(7.1-5)为改进的Euler公式9②用梯形公式：94.几何意义Euler法折线法改进Euler法平均斜率折线法104.几何意义107.1.2截断误差与代数精度定义7.1-1:①称i

=y(xi)–yi为数值解yi的（整体）截断误差。②若yk

=y(xk)，k=0，1，2，…，i–1.由求解公式得数值解，则称为yi

的局部截断误差。注：局部截断误差ei是指单步计算产生的误差，而（整体）截断误差i则考虑到每步误差对下一步的影响。117.1.2截断误差与代数精度11定义7.1-2若求解公式的（整体）截断误差为O(hp)，则称该方法是p阶方法，或是p阶精度。<p阶公式>定理7.1-1设数值解公式：yi+1=yi+h(xi，yi，h)中的函数(x，y，h)关于y满足Lipschitz条件：且其局部截断误差为hp+1阶，则其(整体)截断误差为hp阶，即该数值解公式为p阶方式。注：①局部截断误差较易估计,定理7.1-1表明：若ei

=O(hp+1)则i

=O(hp).②Euler局部截断误差为,所以一阶精度。向后Euler法也是一阶精度。③梯形公式为二阶精度。12定义7.1-2若求解公式的（整体）截断误差为O(hp)，则例1：用Euler方法求解初值问题：取步长h=0.1，并与准确解比较解：因为xi=1+0.1i，而f(x，y)=y+(1+x)y2，故f(xi，yi)=yi+(2+0.1i)yi2于是Euler计算公式为yi+1=yi+0.1[yi+(2+0.1i)yi2]，i=0，1，2，3，413例1：13计算结果如下：

xi1.01.11.21.31.4

y(xi)-1.0000-0.9091-0.8333-0.7692-0.7143

yi

-0.9000-0.8199-0.7540-0.6986-0.6514注：Euler方法精度较低

matlab代码：

yi=-1;y=[yi];fori=1:5y2=yi+0.1*[yi+(2+0.1*i)*yi^2];y=[y,y2];yi=y2;end14计算结果如下：14例2：用改进Euler方法求解初值问题：取步长h=0.1，并与准确解比较解：xi=1+0.1i，于是改进Euler法的计算公式为i=0，1，2，3，4计算结果见P474表7.1-2注：改进Euler方法精度比Euler方法精度高15例2：用改进Euler方法求解初值问题：1516161717xEuler法y改进的Euler法y精确解01.0000001.0000001.0000000.11.0000001.0959091.0954450.21.1918181.1840971.1832160.31.2774381.2662011.2649110.41.3582131.3433601.3416410.51.4351331.4164021.4142140.61.5089661.4859561.48324018xEuler法y改进的Euler法y精确解01.000000xEuler法y改进的Euler法y精确解01.0000001.0000001.0000000.11.0000001.0959091.095445…………0.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.7378671.73205119xEuler法y改进的Euler法y精确解01.0000007.2RungeKutta方法7.2.1构造高阶单步法的直接方法由Taylor公式：当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：207.2RungeKutta方法20当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：即(7.2-1)为p阶方式，上述方式称为Taylor方式。注：利用Tayler公式构造，不实用，高阶导数f

(i)不易计算。21217.2.2RungeKutta方法1.基本思想因为=y(xi)+hf(，y())=y(xi)+hK其中K=f(，y())称为y(x)在[xi，xi+1]上的平均斜率。若取K1=f(xi，y(xi))——Euler公式取K2=f(xi+1，y(xi+1))——向后Euler公式一阶精度取——梯形公式(改进)二阶精度猜想：若能多预测几个点的斜率，再取其加权平均作为K，可望得到较高精度的数值解，从而避免求f的高阶导数。227.2.2RungeKutta方法222.RK公式(7.2-4)其中Kj为y=y(x)在xi+ajh(0

aj

1)处的斜率预测值。aj，bjs，cj为特定常数。232.RK公式233.常数的确定确定的原则是使精度尽可能高。以二阶为例：(7.2-5)(希望y(xi+1)–yi+1=O(hp)的阶数p尽可能高)首先：243.常数的确定24另一方面：由Taylor展开将K2在(xi，yi)处展开。K2=f(xi+a2h，yi+b21hK1

)=f(xi，yi)+a2hf

x'(xi，yi)+b21hK1f

y'(xi，yi)+O(h2).代入(7.2-5)得：yi+1=yi+hc1f(xi，yi)+hc2f(xi，yi)+h2c2[a2f

x'(xi，yi)+b21K1f

y'(xi，yi)]+O(h3)=yi+h(c1+c2)f(xi，yi)+c2a2h2[f

x'(xi，yi)+(b21/a2)f(xi，yi)f

y'(xi，yi)]+O(h3)（希望）25另一方面：由Taylor展开25希望：ei+1=y(xi+1)–yi+1=O(h3).则应：特例：a2=1

c1=c2=1/2，b21=1，得2阶R-K公式：

改进Euler公式。26希望：ei+1=y(xi+1)–yi+1=O(h希望：ei+1=y(xi+1)–yi+1=O(h3).则应：特例：c1=0

c2=1，a2=1/2，b21=1/2，得：

(7.2-7)——中点公式。27希望：ei+1=y(xi+1)–yi+1=O(h4.最常用的R-K公式——标准4阶R-K公式(7.2-8)算法：输入a，b，n，y0h=(b-a)/n，x0=afori=1,i<=n,i++K1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy0=y0+h*(K1+2*K2+2*K3+K4)/6输出x0，y0284.最常用的R-K公式——标准4阶R-K公式输入a，bMatlab代码：functionRunge_Kutta4(a,b,h,y0)n=(b-a)/h;x0=a;fori=1:nK1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy1=y0+h*(K1+2*K2+2*K3+K4)/6;y0=y1end;functionf=f(x,y)f=x+y;end;输入a，b，n，y0h=(b-a)/n，x0=afori=1,i<=n,i++K1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy0=y0+h*(K1+2*K2+2*K3+K4)/6输出x0，y029Matlab代码：输入a，b，n，y0h=(b-a)/n，x例1：（P478）用标准4阶R-K公式求：的数值解。取h=0.2，并与标准解y=2ex

–x–1比较。解：因为f(x，y)=x+y，从而由(7.2-8)得：30例1：（P478）用标准4阶R-K公式求：30例1：（P478）用标准4阶R-K公式求：的数值解。取h=0.2，并与标准解y=2ex

–x–1比较。解：xi+1K1K2K3K4yi+1y(xi+1)0.211.21.221.4441.2428001.2428060.41.4428001.6870801.7115081.9851021.5836361.5836490.61.9836362.2820002.3118362.6460032.0442132.0442380.82.6442133.0086343.0450763.4532282.6510422.6510821.03.4510423.8961463.9406574.4391733.43650331例1：（P478）用标准4阶R-K公式求：xi+1K1K2K注：步长h的选择(P479)使用数值解法求解初值问题，选择步长h是一个重要问题。从每一步看，h小局部截断误差小，整体截断误差也就小；但从整个区间看，h小则节点多，这不仅使计算工作量增大，而且也使舍入误差的累积严重，导致数值不稳定。所以步长h的大小要适当。在实际计算中，常常采用事后估计误差及自动选择步长的办法来保证节点xi处的数值解yi(i=1，2，…，N)满足所要求的精度。32注：步长h的选择(P479)32设所用的一步法是p阶方法，其整体截断误差有渐近公式y(xi+1)–yi+1=Mhp+O(hp+1)其中M是与h无关的常数。考虑用Richardson外推法提高数值解的精度。事实上，从节点xi出发，先以h为步长经一步计算求出数值解，则整体截断误差如上式所示为

(7.2-9)然后以h/2为步长，从节点xi出发经过两步计算求出数值解，则整体截断误差为

(7.2-10)用2p乘(7.2-10)式减去(7.2-9)式，得33设所用的一步法是p阶方法，其整体截断误差有渐近公式33用2p乘(7.2-10)式减去(7.2-9)式，得故有：(7.2-11)(7.2-12)或：(7.2-13)34用2p乘(7.2-10)式减去(7.2-9)式，得34注：(1)(7.2-11)表明若取则整体截断误差为hp+1，即精度提高一阶。(7.2-12)表明数值解的整体截断误差近似为(7.2-13)表明数值解的整体截断误差近似为35注：(1)35(2)记，为精度要求，则当时，有且即步长h/2是合适的；当>

时，说明步长h/2仍然偏大，须将步长减半，继续计算；当<

/2p时，说明已有，步长h/2偏小，应取步长h。36(2)记例2(P481)用变步长的标准4阶R-K方法求初值问题：的数值解，要求精度为37例2(P481)用变步长的标准4阶R-K方法求初值问题：7.3线性多步法希望避免求多个点上f(x，y)的值，并且充分利用前面几步的结果。一般形式：yi+1=a0yi+a1yi-1+…+ap

yi-p+h(b-1y'i+1+b0y'i+…+bp

y'i-p)

(7.3-1)其中y'k

=f(xk，yk)，(k=i–p，i–p+1，…，i，i+1)而ak，bk为待定常数。注：(1)若b–1=0时，(7.3-1)为显式公式，否则为隐式公式。(2)推导方法主要有：数值积分法和Taylor展开法。387.3线性多步法387.3.1数值积分法取节点：xi，xi-1，xi-2，xi-3，作f的三次L-插值多项式有记：fk=f(xk，yk)(k=i，i–1，i–2，i–3)xi-k=xi

–kh，x=xi

+th,(t=0,-1,-2,-3).代入上式有：

397.3.1数值积分法39有

(7.3-2)其中fk=f(xk，yk)(k=i，i–1，i–2，i–3)，称(7.3-2)为4阶Adams显式公式。其局部截断误差(7.3-3)40有40若取节点：xi+1，xi，xi-1，xi-2，作f的3次L-插值多项式，可得4阶Adams隐式公式：

(7.3-4)其局部截断误差：(7.3-5)注：隐式公式的显化：（预测校正）(7.3-6)41若取节点：xi+1，xi，xi-1，xi-2，作f的3次L例：用4阶Adams显式和隐式公式构成的简单预估-校正公式求初值问题：的数值解，取h=0.1。并用标准4阶Runge-Kutta公式求出此初值问题在x1=0.1，x2=0.2，x3=0.3处的数值解，解：计算结果见p487。注：并非所有线性多步法公式(7.3-1)都可用数值积分法得到，但都可用Taylor展开法得到。42例：427.3.2Taylor展开法设yi-k=y(xi

–kh)，y'i-k=y'(xi

–kh)展开为：代入(7.3-1)437.3.2Taylor展开法43代入(7.3-1)得：(7.3-7)44代入(7.3-1)44为使(7.3-1)有m阶精度:y(xi+1)–yi+1=O(hm+1)只须(7.3-7)的前m+1项与y(xi+1)的展式：对应相等，即有方程组:(7.3-8)此时有(7.3-9)45为使(7.3-1)有m阶精度:y(xi+1)–yi+1在(7.3-8)中特别取p=3，m=4有：令a0=a1=a2=b-1=0，可得a3=1，b0=8/3，b1=-4/3，b2=8/3，b3=0.46在46可得a3=1，b0=8/3，b1=-4/3，b2=8/3，b3=0.代入(7.3-1)yi+1=a0yi+a1yi-1+…+ap

yi-p+h(b-1y'i+1+b0y'i+…+bp

y'i-p)

得Milne公式：

(7.3-10)即此时(7.3-11)为4阶精度。47可得a3=1，b0=8/3，b1=-4/3，b2令a1=a3=b2=b3=0，得Hamming公式：(7.3-12)其中：(7.3-13)注：①与单步法相比，多步法不须反复计算f在(xi，xi+1)上某些点处的值，工作量大大减少。②多步法的前几步须用同阶单步法求之。48令a1=a3=b2=b3=0，得Hammin例1:用4阶Adams显式公式求初值问题：的数值解，取h=0.05。解：先用标准4阶Runge-Kutta公式求出此初值问题在x1=0.05，x2=0.1，x3=0.15处的数值解，然后用公式(7.3-2)求其余节点的数值解（p487）49例1:49③Milne显：Hamming隐：显化得Milne-Hamming公式：(7.3-14)50③Milne显：507.4预估校正系统直接预估校正格式：先用显式公式算出预估值，再用同阶隐式公式进行校正，没有充分利用局部截断误差的信息。利用误差补偿的办法，对预估值和校正值进行修正，可以使计算结果的精度更高一些。以4阶Adams为例：显式：局部截断误差(7.3-3)隐式：局部截断误差(7.3-5)517.4预估校正系统51设pi+1，ci+1分别表示xi+1处数值解的预估值和校正值，则(7.4-1)(7.4-2)两式相减：(7.4-3)代入(7.4-1)得(7.4-4)即以作y(xi+1)的预估值，精度提高一阶。52设pi+1，ci+1分别表示xi+1处数值解的预估值和校正值同理，将(7.4-3)代入(7.4-2)得（7.4-5）即以作为y(xi+1)的校正值，其精度可提高一阶。53同理，将(7.4-3)代入(7.4-2)得53注：在(7.4-4)中，由于预估值ci+1还未算出，所以用前一步的ci，pi对pi+1进行修正，由此得：预估修正校正修正(7.4-6)54注：在(7.4-4)中，由于预估值ci+1还未算出，所以用前同理，对Milne公式（显）和Hamming公式（隐）可得带有误差补偿的预估校正系统公式：(7.4-7)注：在(7.4-6)和(7.4-7)中显然无c3，p3，可取p3=c3=0.55同理，对Milne公式（显）和Hamming公式（隐）可得带例1：取步长h=0.1，用带误差补偿的预估——校正公式(7.4-6)求解解：所需y1，y2，y3用标准4阶R-K公式计算，然后按(7.4-6)计算y4，y5.P49256例1：567.5边值问题的差分法基本思想：运用数值微分将导数用离散点上函数值表示，从而将边值问题的微分方程和边界条件转化为只含有限个未知数的差分方程组，并将此差分方程组的解作为该边值问题的数值解。1.二阶常微分方程的第一边值问题(7.5-1)其中q(x)(0)，f(x)在[a，b)上连续，，为常数。577.5边值问题的差分法57设等距节点：xi=a+ih，i=0，1，2，…，n，对其中内节点应用三点微分公式：(7.5-2)i=1，2，…，n–1(当h充分小时，略去O(h2)，并以yi-1，yi，yi+1，代y(xi-1)，y(xi)，y(xi+1)，得计算yi的差分方程组)i=1，2，…，n–1(7.5-3)58设等距节点：xi=a+ih，i=0，1，2，…，加上边界条件即得边值问题(7.5-1)的差分方程组(7.5-4)其中qi

=q(xi)，fi=f(xi).矩阵形式：

(7.5-5)59加上边界条件即得边值问题(7.5-1)的差分方程组59注：①可以证明(7.5-5)存在唯一解，且有误差估计(7.5-6)其中②(7.5-5)等价于：(7.5-5’)为对角占优的三对角方程组，可用追赶法求解。60注：①可以证明(7.5-5)存在唯一解，且有误差估计60例1取h=0.1，求边值问题：的数值解。解：本例的节点xi=0.1i，i=0，1，2，…，10，q(x)=1，f(x)=–x，由(7.5-4)，得差分方程为61例1取h=0.1，求边值问题：61差分方程为由(7.5-5’)得62差分方程为62附：二阶常微分方程三类边界条件设二阶线性常微分方程为常见边界条件有三类：63附：二阶常微分方程三类边界条件设二阶线性常微分方程为63附1差分方程的建立64附1差分方程的建立64656566666767二阶常微分方程边值问题的差分算法68二阶常微分方程边值问题的差分算法686969例题用差分法求解边值问题70例题用差分法求解边值问题707171xyy(x)e=y(x)-y1.00.50.501.10.727985690.72637532-1.610*10-31.21.01403901.0113430-2.696*10-31.31.36782411.3644456-3.378*10-3……………………1.73.68962363.6865656-3.058*10-31.84.56353164.5612288-2.303*10-31.95.58542695.5841425-1.284*10-32.06.77258876.7725887072xyy(x)e=y(x)-y1.00.50.501.10.7数值计算方法第七章常微分方程数值解法73数值计算方法第七章常微分方程数值解法1在工程和科学计算中，所建立的各种常微分方程的初值或边值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因此，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究常微分方程的数值解法是有重要意义的。74在工程和科学计算中，所建立的各种常微分方程的7.0基本概念1.一阶常微分方程的初值问题(7.0-1)注：若f在D={a

x

b,|y|<+}内连续，且满足Lip条件：L

0，使

|f(x，y1)–f(x，y2)|

L|y1–y2|(7.0-2)则(7.0-1)的连续可微解y(x)在[a，b]上唯一存在。757.0基本概念32.初值问题的数值解称(7.0-1)的解y(x)在节点xi处的近似值yi

y(xi)a<x1<x2<...<xn=b.为其数值解，方法称为数值方法。注：①考虑等距节点:xi=a+ih，h=(b–a)/n.②从初始条件y(a)=y0出发，依次逐个计算y1，y2，…，yn的值，称为步进法。两种：单步法、多步法。762.初值问题的数值解4③二阶常微分方程y''(x)=f(x，y(x)，y'

(x))可设为一阶常微分方程组的初值问题：引进新的未知函数z(x)=y'(x)，则其初始条件为：称为一阶微分方程组的初值问题，方法类似。④边界问题，常用差分方法解。77③二阶常微分方程y''(x)=f(x，y(x)，y7.1初值问题数值解法的构造及其精度7.1.1构造方法对于(7.0-1)可借助Taylor展开（导数法）、差商法、积分法实现离散化来构造差分公式：1.Taylor展开：设y

C[a，b]将y(xi+1)=y(xi+h)在xi处展开

[xi，xi+1]

y(xi+1)

yi+hf(xi，yi)其中yi

y(xi).称yi+1=yi+hf(xi，yi).i=0，1，2，...，n–1(7.1-1)为Euler求解公式，（Euler法）787.1初值问题数值解法的构造及其精度62.用差商来表示：得差分方程：

yi+1=yi

+hf(xi，yi).即为Euler公式。若记

yi+1=yi

+hf(xi+1，yi+1).(7.1-2)称为向后Euler法。注：①Euler法为显式，向后Euler法为隐式——须解出yi+1.②可用迭代法yi+1(k+1)=yi

+hf(xi+1，yi+1(k))

k=0，1，2，…解得yi+1，其中yi+1(0)=yi+hf(xi，yi).792.用差商来表示：73.积分法:对(7.0-1)两边取积分得(7.1-3)取不同的数值积分可得不同的求解公式，如：①用矩形公式：左矩形

y(xi+1)

y(xi)+hf(xi，y(xi))Euler公式右矩形

y(xi+1)

y(xi)+hf(xi+1，y(xi+1))向后Euler公式803.积分法:对(7.0-1)两边取积分得8②用梯形公式：

(7.1-4)称(7.1-4)为梯形公式隐式公式。显化：预估值：校正值：(7.1-5)称(7.1-5)为改进的Euler公式81②用梯形公式：94.几何意义Euler法折线法改进Euler法平均斜率折线法824.几何意义107.1.2截断误差与代数精度定义7.1-1:①称i

=y(xi)–yi为数值解yi的（整体）截断误差。②若yk

=y(xk)，k=0，1，2，…，i–1.由求解公式得数值解，则称为yi

的局部截断误差。注：局部截断误差ei是指单步计算产生的误差，而（整体）截断误差i则考虑到每步误差对下一步的影响。837.1.2截断误差与代数精度11定义7.1-2若求解公式的（整体）截断误差为O(hp)，则称该方法是p阶方法，或是p阶精度。<p阶公式>定理7.1-1设数值解公式：yi+1=yi+h(xi，yi，h)中的函数(x，y，h)关于y满足Lipschitz条件：且其局部截断误差为hp+1阶，则其(整体)截断误差为hp阶，即该数值解公式为p阶方式。注：①局部截断误差较易估计,定理7.1-1表明：若ei

=O(hp+1)则i

=O(hp).②Euler局部截断误差为,所以一阶精度。向后Euler法也是一阶精度。③梯形公式为二阶精度。84定义7.1-2若求解公式的（整体）截断误差为O(hp)，则例1：用Euler方法求解初值问题：取步长h=0.1，并与准确解比较解：因为xi=1+0.1i，而f(x，y)=y+(1+x)y2，故f(xi，yi)=yi+(2+0.1i)yi2于是Euler计算公式为yi+1=yi+0.1[yi+(2+0.1i)yi2]，i=0，1，2，3，485例1：13计算结果如下：

xi1.01.11.21.31.4

y(xi)-1.0000-0.9091-0.8333-0.7692-0.7143

yi

-0.9000-0.8199-0.7540-0.6986-0.6514注：Euler方法精度较低

matlab代码：

yi=-1;y=[yi];fori=1:5y2=yi+0.1*[yi+(2+0.1*i)*yi^2];y=[y,y2];yi=y2;end86计算结果如下：14例2：用改进Euler方法求解初值问题：取步长h=0.1，并与准确解比较解：xi=1+0.1i，于是改进Euler法的计算公式为i=0，1，2，3，4计算结果见P474表7.1-2注：改进Euler方法精度比Euler方法精度高87例2：用改进Euler方法求解初值问题：1588168917xEuler法y改进的Euler法y精确解01.0000001.0000001.0000000.11.0000001.0959091.0954450.21.1918181.1840971.1832160.31.2774381.2662011.2649110.41.3582131.3433601.3416410.51.4351331.4164021.4142140.61.5089661.4859561.48324090xEuler法y改进的Euler法y精确解01.000000xEuler法y改进的Euler法y精确解01.0000001.0000001.0000000.11.0000001.0959091.095445…………0.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.7378671.73205191xEuler法y改进的Euler法y精确解01.0000007.2RungeKutta方法7.2.1构造高阶单步法的直接方法由Taylor公式：当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：927.2RungeKutta方法20当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：即(7.2-1)为p阶方式，上述方式称为Taylor方式。注：利用Tayler公式构造，不实用，高阶导数f

(i)不易计算。93217.2.2RungeKutta方法1.基本思想因为=y(xi)+hf(，y())=y(xi)+hK其中K=f(，y())称为y(x)在[xi，xi+1]上的平均斜率。若取K1=f(xi，y(xi))——Euler公式取K2=f(xi+1，y(xi+1))——向后Euler公式一阶精度取——梯形公式(改进)二阶精度猜想：若能多预测几个点的斜率，再取其加权平均作为K，可望得到较高精度的数值解，从而避免求f的高阶导数。947.2.2RungeKutta方法222.RK公式(7.2-4)其中Kj为y=y(x)在xi+ajh(0

aj

1)处的斜率预测值。aj，bjs，cj为特定常数。952.RK公式233.常数的确定确定的原则是使精度尽可能高。以二阶为例：(7.2-5)(希望y(xi+1)–yi+1=O(hp)的阶数p尽可能高)首先：963.常数的确定24另一方面：由Taylor展开将K2在(xi，yi)处展开。K2=f(xi+a2h，yi+b21hK1

)=f(xi，yi)+a2hf

x'(xi，yi)+b21hK1f

y'(xi，yi)+O(h2).代入(7.2-5)得：yi+1=yi+hc1f(xi，yi)+hc2f(xi，yi)+h2c2[a2f

x'(xi，yi)+b21K1f

y'(xi，yi)]+O(h3)=yi+h(c1+c2)f(xi，yi)+c2a2h2[f

x'(xi，yi)+(b21/a2)f(xi，yi)f

y'(xi，yi)]+O(h3)（希望）97另一方面：由Taylor展开25希望：ei+1=y(xi+1)–yi+1=O(h3).则应：特例：a2=1

c1=c2=1/2，b21=1，得2阶R-K公式：

改进Euler公式。98希望：ei+1=y(xi+1)–yi+1=O(h希望：ei+1=y(xi+1)–yi+1=O(h3).则应：特例：c1=0

c2=1，a2=1/2，b21=1/2，得：

(7.2-7)——中点公式。99希望：ei+1=y(xi+1)–yi+1=O(h4.最常用的R-K公式——标准4阶R-K公式(7.2-8)算法：输入a，b，n，y0h=(b-a)/n，x0=afori=1,i<=n,i++K1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy0=y0+h*(K1+2*K2+2*K3+K4)/6输出x0，y01004.最常用的R-K公式——标准4阶R-K公式输入a，bMatlab代码：functionRunge_Kutta4(a,b,h,y0)n=(b-a)/h;x0=a;fori=1:nK1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy1=y0+h*(K1+2*K2+2*K3+K4)/6;y0=y1end;functionf=f(x,y)f=x+y;end;输入a，b，n，y0h=(b-a)/n，x0=afori=1,i<=n,i++K1=f(x0,y0)K2=f(x0+h/2,y0+h*K1/2)K3=f(x0+h/2,y0+h*K2/2)K4=f(x0+h,y0+h*K3)x0=x0+hy0=y0+h*(K1+2*K2+2*K3+K4)/6输出x0，y0101Matlab代码：输入a，b，n，y0h=(b-a)/n，x例1：（P478）用标准4阶R-K公式求：的数值解。取h=0.2，并与标准解y=2ex

–x–1比较。解：因为f(x，y)=x+y，从而由(7.2-8)得：102例1：（P478）用标准4阶R-K公式求：30例1：（P478）用标准4阶R-K公式求：的数值解。取h=0.2，并与标准解y=2ex

–x–1比较。解：xi+1K1K2K3K4yi+1y(xi+1)0.211.21.221.4441.2428001.2428060.41.4428001.6870801.7115081.9851021.5836361.5836490.61.9836362.2820002.3118362.6460032.0442132.0442380.82.6442133.0086343.0450763.4532282.6510422.6510821.03.4510423.8961463.9406574.4391733.436503103例1：（P478）用标准4阶R-K公式求：xi+1K1K2K注：步长h的选择(P479)使用数值解法求解初值问题，选择步长h是一个重要问题。从每一步看，h小局部截断误差小，整体截断误差也就小；但从整个区间看，h小则节点多，这不仅使计算工作量增大，而且也使舍入误差的累积严重，导致数值不稳定。所以步长h的大小要适当。在实际计算中，常常采用事后估计误差及自动选择步长的办法来保证节点xi处的数值解yi(i=1，2，…，N)满足所要求的精度。104注：步长h的选择(P479)32设所用的一步法是p阶方法，其整体截断误差有渐近公式y(xi+1)–yi+1=Mhp+O(hp+1)其中M是与h无关的常数。考虑用Richardson外推法提高数值解的精度。事实上，从节点xi出发，先以h为步长经一步计算求出数值解，则整体截断误差如上式所示为

(7.2-9)然后以h/2为步长，从节点xi出发经过两步计算求出数值解，则整体截断误差为

(7.2-10)用2p乘(7.2-10)式减去(7.2-9)式，得105设所用的一步法是p阶方法，其整体截断误差有渐近公式33用2p乘(7.2-10)式减去(7.2-9)式，得故有：(7.2-11)(7.2-12)或：(7.2-13)106用2p乘(7.2-10)式减去(7.2-9)式，得34注：(1)(7.2-11)表明若取则整体截断误差为hp+1，即精度提高一阶。(7.2-12)表明数值解的整体截断误差近似为(7.2-13)表明数值解的整体截断误差近似为107注：(1)35(2)记，为精度要求，则当时，有且即步长h/2是合适的；当>

时，说明步长h/2仍然偏大，须将步长减半，继续计算；当<

/2p时，说明已有，步长h/2偏小，应取步长h。108(2)记例2(P481)用变步长的标准4阶R-K方法求初值问题：的数值解，要求精度为109例2(P481)用变步长的标准4阶R-K方法求初值问题：7.3线性多步法希望避免求多个点上f(x，y)的值，并且充分利用前面几步的结果。一般形式：yi+1=a0yi+a1yi-1+…+ap

yi-p+h(b-1y'i+1+b0y'i+…+bp

y'i-p)

(7.3-1)其中y'k

=f(xk，yk)，(k=i–p，i–p+1，…，i，i+1)而ak，bk为待定常数。注：(1)若b–1=0时，(7.3-1)为显式公式，否则为隐式公式。(2)推导方法主要有：数值积分法和Taylor展开法。1107.3线性多步法387.3.1数值积分法取节点：xi，xi-1，xi-2，xi-3，作f的三次L-插值多项式有记：fk=f(xk，yk)(k=i，i–1，i–2，i–3)xi-k=xi

–kh，x=xi

+th,(t=0,-1,-2,-3).代入上式有：

1117.3.1数值积分法39有

(7.3-2)其中fk=f(xk，yk)(k=i，i–1，i–2，i–3)，称(7.3-2)为4阶Adams显式公式。其局部截断误差(7.3-3)112有40若取节点：xi+1，xi，xi-1，xi-2，作f的3次L-插值多项式，可得4阶Adams隐式公式：

(7.3-4)其局部截断误差：(7.3-5)注：隐式公式的显化：（预测校正）(7.3-6)113若取节点：xi+1，xi，xi-1，xi-2，作f的3次L例：用4阶Adams显式和隐式公式构成的简单预估-校正公式求初值问题：的数值解，取h=0.1。并用标准4阶Runge-Kutta公式求出此初值问题在x1=0.1，x2=0.2，x3=0.3处的数值解，解：计算结果见p487。注：并非所有线性多步法公式(7.3-1)都可用数值积分法得到，但都可用Taylor展开法得到。114例：427.3.2Taylor展开法设yi-k=y(xi

–kh)，y'i-k=y'(xi

–kh)展开为：代入(7.3-1)1157.3.2Taylor展开法43代入(7.3-1)得：(7.3-7)116代入(7.3-1)44为使(7.3-1)有m阶精度:y(xi+1)–yi+1=O(hm+1)只须(7.3-7)的前m+1项与y(xi+1)的展式：对应相等，即有方程组:(7.3-8)此时有(7.3-9)117为使(7.3-1)有m阶精度:y(xi+1)–yi+1在(7.3-8)中特别取p=3，m=4有：令a0=a1=a2=b-1=0，可得a3=1，b0=8/3，b1=-4/3，b2=8/3，b3=0.118在46可得a3=1，b0=8/3，b1=-4/3，b2=8/3，b3=0.代入(7.3-1)yi+1=a0yi+a1yi-1+…+ap

yi-p+h(b-1y'i+1+b0y'i+…+bp

y'i-p)

得Milne公式：

(7.3-10)即此时(7.3-11)为4阶精度。119可得a3=1，b0=8/3，b1=-4/3，b2令a1=a3=b2=b3=0，得Hamming公式：(7.3-12)其中：(7.3-13)注：①与单步法相比，多步法不须反复计算f在(xi，xi+1)上某些点处的值，工作量大大减少。②多步法的前几步须用同阶单步法求之。120令a1=a3=b2=b3=0，得Hammin例1:用4阶Adams显式公式求初值问题：的数值解，取h=0.05。解：先用标准4阶Runge-Kutta公式求出此初值问题在x1=0.05，x2=0.