1、第五章 常微分方程数值解法 (Numerical Methods for Ordinary Differential Equations),引言,引言, 考虑一阶常微分方程的初值问题 (Initial-Value Problem ):,只要 f (x, y) 在x a, b 上连续，且关于 y 满足 Lipschitz 条件，即存在与 x, y 无关的常数 L 使对任意定义在 a, b 上的 y1(x) 和 y2(x) 都成立，则上述常微分方程存在唯一解。,引言,数值解法就是要计算出解函数 y(x) 在一系列节点 a = x0 x1 xn= b 处的近似值y0 y1 yn,节点间距 称为步长，
2、通常采用等距节点，即取 hi = h (常数)。,解决： 数值解法的一个基本特点是“步进式”，即求解时顺着节点排列的次序一步步地向前推进。,注意：与“迭代法”区别,第五章 常微分方程数值解法,欧拉方法（ Eulers Method ）,显式Euler公式,由两点公式求导数，在xj , xj+1子区间上有：,其中j xj , xj+1,代入方程 有,显式Euler公式,显式Euler公式的误差,局部截断误差,显式Euler公式,（P169）例1：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题,依次下去计算结果见P169,隐式Euler公式,由两点公式求
3、导数，在xj , xj+1子区间上有：,其中j xj , xj+1,代入方程 有,局部截断误差,隐式Euler公式,是一个关于yj+1的方程, 要从中解出yj+1,（P169）例1：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题,计算过程见P170,用隐式Euler法求解，有：,解：,从中解出yi+1，有：,显隐结合的预测校正系统避免求解方程 （predictor-corrector method ）,Step 1: 先用显式欧拉公式作预测，算出预测值,Step 2: 再用隐式欧拉公式作校正，得到校正值,写成一个公式为：,计算顺序：,（P170）例1
4、：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题,计算过程见P170,解：,用预测校正系统求解，有：,y0=1,p阶精度, 欧拉法的局部截断误差：,欧拉法具有 1 阶精度。,显式：,隐式：,梯形公式 （trapezoid formula）,注：的确有局部截断误差 ，即梯形公式具有2 阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式。,从显式、隐式Euler法和Euler局部截断误差来看似乎可以有如下式子梯形公式：,隐式Euler,显式Euler,基于数值积分的求解思想,求出 积分的近似表达， 也就求出了y(xk),h,xk,x,f(x,y),f(
5、x,y),xk-1,f(xk,yk),f(xk-1,yk-1),基于数值积分的求解思想,用矩形求积公式代替,左矩形显式Euler公式,右矩形隐式Euler公式,基于数值积分的求解思想,用梯形求积公式代替,注：梯形公式是隐式公式显式、隐式Euler的平均梯形公式具有2阶精度,余项,Cant you give me a formula with all the advantages yet without any of the disadvantages?,Let me try！,Euler公式及梯形公式总结,改进的Euler法预测校正系统 （ modified Eulers method）,注：
6、可以证明该算法具有 2 阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。它的稳定性高于显式欧拉法。误差比较：梯形法改进Euler法 显、隐Euler法,Step 1: 先用显式欧拉公式作预测，算出预测值,Step 2: 再用梯形公式作校正，得到校正值,举例,(P173例3)：取h=0.1，分别用Euler法、梯形公式、改进的Euler法求解初值问题,解：,计算过程参见P173,Euler法,梯形法,解出yk得：,举例,(P173例3)：取h=0.1，分别用Euler法、梯形公式、改进的Euler法求解初值问题,解：,计算过程参见P173,改进Euler法,Simpson公
7、式（Simpson formula ）,注：Simpson公式是一个隐式公式Simpson公式是多步法：yk-1, yk-2=yk,应用simpson求积公式：,Simpson公式,4阶精度,P198 习题6 用Euler法和改进的Euler法解方程，h0.1 注意改题目为0=x=0.4,第五章 作 业,欧拉方法（ Eulers Method ）,第五章 常微分方程数值解法,龙格库塔法 （RungeKutta Method ）,二阶RungeKutta法改进的Euler法, 考察改进的欧拉法，可以将其改写为：,K1,K2,xi+h,步长一定是一个h 吗？,斜率一定 取K1 K2 的 平均值吗？
8、,龙格 - 库塔法（ RungeKutta Method ）,首先希望能确定系数 1、2、p，使得到的算法格式有2阶精度，即在 的前提假设下，使得,Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开,龙格 - 库塔法（ RungeKutta Method ）,Step 2: 将 K1 K2 代入第1式，得到,Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较,要求 ，则必须有：,龙格 - 库塔法（ RungeKutta Method ）,这里有 个未知数， 个方程。,3,2,存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格
9、式。,注意: 就是改进的欧拉法, 也就是2阶R-K法。, : 为获得更高的精度，应该如何进一步推广？,三阶龙格 - 库塔法,三、四阶R-K法思想：利用各种显化方法，对以上公式显化,xi,xi +h,xi +h/2,三阶龙格 - 库塔法,K1,K2,K3,显Euler,P177,要计算3个函数 可以证明(Taylor展开K1，K2，K3）e(h)=O(h4)具有3阶精度,P180 例：取h=0.1，用三阶、四阶R-K法求解初值问题,举例,解：用三阶R-K法,取n=1 计算K1， K2， K3， = y1,取n=2 计算K1， K2， K3， = y2,取n=3 计算K1， K2， K3， = y
10、3,列表计算,四阶龙格 - 库塔法(经典R-K法) /* Classical Runge-Kutta Method */,e(h)=O(h5) 4阶精度,优缺点： 优点是：(1)都是一步法，因此只要给定一个初始值就可以一直计算下去；(2)精度相对较高，如经典R-K法为四阶精度 缺点是：(特别是三阶和四阶法)计算量较大。,其他问题讨论： R-K法推导基于Taylor展式，因而要求y(x)有较好的光滑性（即有高阶导数）。 最常用的是四阶公式，它适用于一般的问题，准确、稳定、易于编程。 步长h减小，局部误差O(h5)减小但步数增加，舍入误差积累增加,四阶龙格 - 库塔法(经典R-K法) /* Cla
11、ssical Runge-Kutta Method */,P180 例：取h=0.1，用三阶、四阶R-K法求解初值问题,解：用三阶R-K法,取n=1 计算K1， K2， K3， K4 = y1,取n=2 计算K1， K2， K3， K4 = y2,取n=3 计算K1， K2， K3， K4 = y3,列表计算,例：取h=0.2，用四阶龙格库塔法求解初值问题,举例,解：,这里 ，经典的四阶龙格-库塔公式为,举例,表中列出了计算结果，同时列出了相应的精确解.比较本章第一个例子的计算结果，显然龙格-库塔方法的精度高.,名称RK4（A,B,Y0,H） 存储：K1K4 = K1K4 YN = y1,y2
12、yn 算法：(1) 输入a，b，n，y0，设计函数F(x,y)(2)Y0=y0, H=(b-a)/n，H1=H/2(3)对i=0，1，2， ,N-1(由Yi计算Yi+1)xi=x0+i*hK1=F(xi,Yi)K2=F(xi+H1,Yi+H1*K1)K3=F(xi+H1,Yi+H1*K2)K4=F(xi+H,Yi+H*K3)Yi+1=Yi+(K1+2*K2+2*K3+K4)*h/6(4)输出Y1, Y2,YN.,四阶龙格 - 库塔法算法设计（P182),P199 习题10（2）,第五章 作 业,龙格库塔法 （RungeKutta Method ）,第五章 常微分方程数值解法,线性多步法（ Mu
13、ltistep Method ）,线性多步法（ Multistep Method ）,对应xk： a = x0 x1 xn= b的yk： y0 y1 yn,优点：不需要其他计算值，计算yk只要用到yk-1 缺点：没有充分利用前几步得到的信息 处理y(xk-1+ah)的困难,用到p步计算的结果,开型求解Adams法,思想：计算积分时，插值节点取在积分区间xk-1,xk的外部 三阶Adams： 显式：取xk-1, xk-2, xk-3 作节点f(x,y)L2(x) 误差O(h4) 隐式：取xk, xk-1, xk-2 作节点f(x,y)L2(x) 误差O(h4) 四阶Adams： 显式：取xk-1
14、, xk-2, xk-3 , xk-4 作节点f(x,y)L3(x) 误差O(h5) 隐式：取xk , xk-1, xk-2, xk-3 作节点f(x,y)L3(x) 误差O(h5),预测校正系统：显式Adams隐式Adams（三阶、四阶）,三阶显式Adams法,取xk-1, xk-2, xk-3 作节点f(x,y)L2(x) 误差O(h4),简记为：,误差：,3阶精度,三阶隐式Adams法,取xk, xk-1, xk-2 作节点f(x,y)L2(x) 误差O(h4),误差：,3阶精度,四阶显式Adams法,取xk-1, xk-2, xk-3 , xk-4 作节点f(x,y)L3(x) 误差O
15、(h5),误差：,4阶精度,四阶隐式Adams法,取xk , xk-1, xk-2, xk-3 作节点f(x,y)L3(x) 误差O(h5),误差：,4阶精度,亚当姆斯预测-校正系统 （ Adams predictor-corrector system ）,四阶误差：,三阶,四阶,其中：,亚当姆斯预测-校正系统 （ Adams predictor-corrector system ）,Step 1: 用Runge-Kutta 法计算前 k 个初值；,Step 2: 用Adams 显式计算预测值；,Step 3: 用同阶Adams 隐式计算校正值。,注意：三步所用公式的精度必须相同。通常用经典R
16、unge-Kutta 法（4阶）配合4阶Adams 公式。,闭型求解Milne公式,闭型求解思想：,扩大积分区间：xk-1,xk xk-4,xk,增加内部节点： xk-1, xk-2, xk-3 , xk-4 原来的外部节点 积分区间的内部节点,方程形式变化：,闭型求解Milne公式,推导： (1) 取xk-1, xk-2, xk-3 , xk-4 为节点作三次Newton前插公式代替f(x, y(x)：f(x, y(x)= N3(x)+R4(x) (2) xxk-4+th 则 xxk-4,xk时，t 0,4 作x-t替换 (3) 计算积分及误差,截断误差：,4阶精度,Milne公式,闭型求解
17、公式,Milne公式不能直接使用，一般是首先用四阶RK法求出前几个节点的值，然后使用Milne公式求解。 前面知道Simpson公式（四阶精度）是隐式积分公式，所以利用Milne公式，可以得到如下的精度为四阶的预测校正系统。,解： (1) 用三阶R-K法求y(0.4) 的值,P187例：分别用三阶Adams显式公式和Milne公式计算下面的初值 问题的解y(x)在0.6和0.8处的值， h0.2，已知y1=0.181。,解： (1) 用三阶R-K法求y(0.4) 的值,P187例：分别用三阶Adams显式公式和Milne公式计算下面的初值 问题的解y(x)在0.6和0.8处的值， h0.2，已
18、知y1=0.181。,(2) 将0、0.2、0.4处的值代入三阶显式Adams公式计算y(0.6),(3) 将0、0.2、0.4、0.6处的值代入Milne公式计算y(0.8),计算过程见P187,第五章 常微分方程数值解法,线性多步法（ Multistep Method ）,P198 习题15（第三版） P200 习题16（第二版）,第五章 常微分方程数值解法,欧拉方法（ Eulers Method ）,例:用改进的Euler法推出Yn和Xn的关系表达式，并与准确解析式作比较,举例,补充：,导出改进Euler公式的解表达并与,解：,其中 为Euler算出,Euler公式：,第五章 常微分方程数值解法,龙格库塔法 （RungeKutta Method ）,例:用四阶R-K法求解初值问题,举例,解：,举例,举例,举例,例:用四阶Adams求解初值问题,举例,解：,举例,第五章 常微分方程数值解法,线性多步法（ Multistep Method ）,
