计算方法常微分方程的数值解法

1、2020/7/6,1,第九章 常微分方程的数值解法,9.1、引言 9.2、初值问题的数值解法-单步法 9.3、龙格-库塔方法 9.4、收敛性与稳定性,主 要 内 容,2020/7/6,2,主 要 内 容,研究的问题 数值解法的意义,9.1 引 言,2020/7/6,3,现实世界中大多数事物,内部联系非常复杂,找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系,此种关系的数学表达就为,微分方程,1.什么是微分方程 ?,其状态随着 时间、地点、条件 的不同而不同,2020/7/6,4,如何利用数值方法求解微分方程（组）的问题。,2.数值求解微分方程的意义,如何建立数学

2、模型？ 各类微分方程本身和他们的解所具有的特性 已在常微分方程及数学物理方程中得以解释，,本章专门讨论,2020/7/6,5,常微分方程表示方法,2020/7/6,6,本章主要讨论一阶常微方程的初值问题,基本知识,其中f (x,y)是已知函数，(1.2)是定解条件也称为 初值条件。,各种数值解法,2020/7/6,7,则称 f (x,y) 对y 满足李普希兹条件，L 称为Lipschitz常数.,常微分方程的理论指出：,当 f (x,y) 定义在区域 G=(axb,y),若存在正的常数 L 使：,就可保证方程解的存在唯一性,(Lipschitz)条件,2020/7/6,8,可以证明,如果函数在

3、带形区域 R=axb, -y内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使,对R内任意两个 都成立,则方程( 7.1 )的解 在a, b上存在且唯一。,2020/7/6,9,寻找解析解的过程称为求解微分方程组。,一个或一组具有所要求阶连续导数的解析函数，将它代入微分方程（组），恰使其所有条件都得到满足的解称为解析解(或古典解),称为真解或解。,3.什么是微分方程 (组)的解析解?,在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。,2020/7/6,10,当

4、x=0时,y=1,可得c=1特解,当x=0时,y=1,可得c=-1特解,两边积分,通解,常微分方程解法回顾,-,2020/7/6,11,但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。 譬如,这个一阶微分方程就不能用初等函数及其积分来表达它的解。,3.什么是微分方程 (组)的解析解?,再如，方程,的解 ,虽然有表可查,但对于表上没有给出 的值,仍需插值方法来计算,2020/7/6,12,在微分方程中, 自变量的个数只有一个, 称为常微分方程. 自变量的个数为两个或两个以上的微分方程叫偏微分方程。,2020/7/6,13,4.什么是微分方程的数值解?,虽然求解微分方程有许

5、多解析方法，但解析方法只能够求解一些特殊类型的方程，从实际意义上来讲，我们更关心的是某些 特定的自变量在某一个定义范围内的一系列离散点上的近似值。,寻找数值解的过程称为数值求解微分方程。,把这样一组近似解称为 微分方程在该范围内的,数值解,2020/7/6,14,在大量的实际方程中出现的函数起码的连续性都无法保证，更何况要求阶的导数,求解数值解,很多微分方程 根本求不到 问题的解析解！,重要手段。,2020/7/6,15,常微分方程的数值解法常用来求近似解,根据提供的算法,通过计算机,便捷地实现,5.常微分方程数值解法的特点,数值解法得到的近似 解(含误差)是一个 离散的函数表.,2020/7

6、/6,16,我们以下的讨论，都在满足上述条件下进行。,若 f (x,y) 在区域 G连续，关于y,一阶常微分方程的初值问题的解存在且唯一。,满足李普希兹 条件,一阶常微分方程组常表述为：,方程组,初值条件,2020/7/6,17,写成向量形式为,高阶常微分方程定解问题如二阶定解问题：,2020/7/6,18,也就解决了高阶方程的定解问题。,这些解法都可以写成向量形式,用于一阶常微分方程组的初值问题。,2020/7/6,19,简单的数值方法与基本概念,1. 简单欧拉法（Euler） 2后退的欧拉法 3梯形法 4改进Euler法,9.2、初值问题的数值解法单步法,2020/7/6,20,1. 简单

7、的欧拉(Euler)方法,考虑模型：,在精度要求不高时使用。,通过欧拉方法的讨论,弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路.,欧拉方法,最简单而直观 实用方法,2020/7/6,21,2. 欧拉方法的导出,把区间a,b,分为n个小区间,步长为,要计算出解函数 y(x) 在一系列节点,节点,处的近似值,N等分,2020/7/6,22,对微分方程(1.1)两端从,进行积分,2020/7/6,23,右端积分用左矩形数值求积公式,得,或用向前差 商近似导数,2020/7/6,24,亦称为欧拉折线法 /* Eulers polygonal arc method*/,依上述公式逐次计算可

8、得：,故也称Euler为单步法。,公式右端只含有已知项,所以又称为显格式的单步法。,每步计算,只用到,2020/7/6,25,3.欧拉公式有明显的几何意义,依此类推得到一折线,2020/7/6,26,就是用这条折线近似地代替曲线,欧拉方法,从几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙。,2020/7/6,27,利用 Euler方法求初值问题,解 此时的 Euler公式为,例,的数值解.此问题的精确解是,分别取步长h=0.2 ,0.1 ,0.05,计算结果如下,2020/7/6,28,2020/7/6,29,4.欧拉法的局部截断误差：,如果单步差分公式的局部截

9、断误差为O(hp+1)，则称该公式为p阶方法。这里p为非负整数。显然,阶数越高，方法的精度越高。,2020/7/6,30,Ri 的主项 /* leading term */, 欧拉法的局部截断误差：,欧拉法具有 1 阶精度。,Taylor展开式,一元函数的Taylor展开式为:,2020/7/6,31,5. 欧拉公式的改进：, 隐式欧拉法 /* implicit Euler method */,由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。,一般先用显式计算一个初值，再迭代求

10、解。, 隐式欧拉法的局部截断误差：,即隐式欧拉公式具有 1 阶精度。,见P285,2020/7/6,32,6.梯形公式 /* trapezoid formula */, 显、隐式两种算法的平均,注：的确有局部截断误差 ， 即梯形公式具有2 阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。, 中点欧拉公式 /* midpoint formula */,假设 ，则可以导出 即中点公式具有 2 阶精度。,需要2个初值 y0和 y1来启动递推 过程，这样的算法称为双步法 /* double-step method */，而前面的三种算法都是单步

11、法 /* single-step method */。,2020/7/6,33,方 法,显式欧拉,隐式欧拉,梯形公式,中点公式,简单,精度低,稳定性最好,精度低, 计算量大,精度提高,计算量大,精度提高, 显式,多一个初值, 可能影响精度,2020/7/6,34,例 用欧拉方法，隐式欧拉方法和欧拉中点公式计算,的近似解，取步长h=0.1，并与精确值比较,解:欧拉方法的算式为：,欧拉隐式方法在本题可解出方程，不必迭代，公式为：,欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式,本题精确解为y=e-x,2020/7/6,35, 改进欧拉法 /* modified Eulers method

12、*/,注：此法亦称为预测-校正法 /* predictor-corrector method */。可以证明该算法具有 2 阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。,2020/7/6,36,例 用欧拉公式和梯形公式的预估校正法（改进欧拉公式）计算：,的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较.方程的解析解为:,解: 本例中欧拉公式为：,2020/7/6,37,梯形公式只校正一次的格式为,结果列入表9.2,2020/7/6,38,改进欧拉法算法实现 （1）计算步骤 输入 , h , N 使用以下改进欧拉法公式进行

13、计算 输出 ，并使 (2) 转到 直至n N 结束。,2020/7/6,39,改进欧拉法的流程图,2020/7/6,40,9.3 龙格 - 库塔法 /* Runge-Kutta Method */,斜率 一定取K1 K2 的平均值吗？,步长一定是一个h 吗？,单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。,建立高精度的单步递推格式。,2020/7/6,41,首先希望能确定系数 1、2、p，使得到的算法格式有2阶精度，即在 的前提假设下，使得,Step 1: 将 K2 在 ( x

14、i , yi ) 点作 Taylor 展开,Step 2: 将 K2 代入第1式，得到,2020/7/6,42,Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较,要求 ，则必须有：,这里有 个未知数， 个方程。,3,2,存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格式。,注意到， 就是改进的欧拉法。,Q: 为获得更高的精度，应该如何进一步推广？,2020/7/6,43,构造高阶单步法的直接方法 由Taylor公式： 当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程： 其局部截断误差为：,差

15、分方程为p阶方式，上述方式称为Taylor方式。 注：利用Tayler公式构造，不实用，高阶导数f (i)不易计算。,2020/7/6,44,RungeKutta方法 1. 基本思想 因为 其中K=f(,y()称为y(x)在xi,xi+1上的平均斜率。 若取K1=f(xi,y(xi)Euler公式 取K2=f(xi+1,y(xi+1)向后Euler公式 一阶精度 取 梯形公式 二阶精度 猜想：若能多预测几个点的斜率，再取其加权平均作为K，可望得到较高精度的数值解，从而避免求f 的高阶导数。,2020/7/6,45,2. RK公式 其中Kj为y = y(x)在xi + ajh (0 aj 1)处

16、的斜率预测值。 aj，bjs，cj为特定常数。,2020/7/6,46,3. 常数的确定 确定的原则是使精度尽可能高。 以二阶为例： (希望y(xi+1) yi+1 = O(hp)的阶数p尽可能高) 首先：,2020/7/6,47,另一方面： 将K2在(xi，yi)处展开。 K2 = f (xi,yi) + a2hf x(xi,yi) + b21hK1 f y(xi,yi) + O(h2) 可得： yi+1 = yi + hc1 f (xi,yi) + hc2 f (xi,yi) + h2c2a2 f x(xi,yi) + b21K1 f y(xi,yi) + O(h3) = yi + h(c

17、1 + c2) f (xi,yi) + c2a2h2f x(xi,yi) + (b21/a2) f (xi,yi) f y(xi,yi)+O(h3) （希望）,2020/7/6,48,希望：ei+1 = y(xi+1) yi+1 = O(h3). 则应： 特例：a2 = 1 c1 = c2 = 1/2，b21 = 1，得2阶R-K公式： 改进欧拉公式。,2020/7/6,49,希望：ei+1 = y(xi+1) yi+1 = O(h3). 则应： 特例：c1 = 0 c2 = 1，a2 = 1/2，b21 = 1/2，得： 称为中点公式。,2020/7/6,50,4. 最常用的R-K公式 标准

18、4阶R-K公式 算法：,2020/7/6,51,Matlab代码： function Runge_Kutta4(a,b,h,y0) n=(b-a)/h;x0=a; for i=1:n 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+h y1=y0+h*(K1+2*K2+2*K3+K4)/6; y0=y1 end; end; function f=f(x,y) f=x+y; end;,2020/7/6,52,例：用标准4阶R-K公式求： 的数值解。取h = 0.2，并与标准解y

19、= 2ex x 1比较。 解：因为f(x，y) = x + y，从由得：,2020/7/6,53,2020/7/6,54, 由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h 取小。,深入研究龙格- 库塔法请看此处！,2020/7/6,55,9.4 收敛性与稳定性 /* Convergency and Stability */, 1.收敛性 /* Convergency */,例：就初值问题 考察欧拉显式格式的收敛性。,解：该问题的精确解为,欧拉公式为,对任意固定的 x = xi = i h ，有,2020/7/6,56, 2.

20、稳定性 /* Stability */,例：考察初值问题 在区间0, 0.5上的解。 分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。,1.0000 2.0000 4.0000 8.0000 1.6000101 3.2000101,1.0000 2.5000101 6.2500102 1.5625102 3.9063103 9.7656104,1.0000 2.5000 6.2500 1.5626101 3.9063101 9.7656101,1.0000 4.9787102 2.4788103 1.2341104 6.1442106 3.0590107,What is wrong ?!,2

21、020/7/6,57,一般分析时为简单起见，只考虑试验方程 /* test equation */,常数，可以是复数,2020/7/6,58,例：考察隐式欧拉法,可见绝对稳定区域为：,注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。,2020/7/6,59,例：隐式龙格-库塔法,而显式 1 4 阶方法的绝对稳定区域为,其中2阶方法 的绝对稳定区域为,无条件稳定,2020/7/6,60,1. Runge-Kutta 法的一般形式 2. 2阶Runge-Kutta 方法 3. 经典Runge-Kutta 方法 4R-K-Fehhlberg 方法 5. 隐式R-K方法 6. 变步长方法,龙格

22、库塔法深入研究,2020/7/6,61,1.Runge-Kutta 法的一般形式,Runge-Kutta 法是用区间上若干个点上的导数的线性组合得到平均斜率，以提高方法的阶。其一般形式为 ：,2020/7/6,62,式(9.11) 称L级p阶Runge-Kutta方法(简称R-K法)。 当L1就是欧拉法，当L2时为改进的欧拉法。,其中,它的局部截断误差是,(9.11),2. 2级2阶Runge-Kutta方法,令 L=2，则,2020/7/6,63,其局部截断误差是：,这是3个未知量的两个方程，其中有一个自由参数，方程组有无穷多解，从而得到一族2级2阶R-K方法 。,2020/7/6,64,3

23、. 经典Runge-Kutta方法,我们可以构造出一族3级3阶，一族4级4阶和一族5级4 阶等R-K方法。最常用的4级4阶是如下的经典R-K方法：,2020/7/6,65,4R-K-Fehhlberg 方法,R-K-Fehhlberg方法是在R-K方法的基础上引进误差和步长 控制的办法。即利用5阶R-K法,估计4阶R-K的局部误差，其中,2020/7/6,66,注：R-K-Fehhlberg比4阶R-K方法具有更大的优越性，他是计算稳定，高精度的方法，他的显著优点是，每一步仅需计算f的6个值；若用4阶R-K-L与5阶R-K-L在一起使用，则每步需要计算f的10个值。推荐使用！,2020/7/6

24、,67,5. 隐式R-K方法,类似于显式R-K公式，稍加改变，就得到隐式R-K方法。,它与显式R-K公式的区别在于：显式公式中对系数求和的上限是i-1，从而构成的矩阵是一个严格下三角阵。而在隐式公式中对系数求和的上限是L，从而构成的矩阵是方阵，需要用迭代法求出Ki，。推导隐式公式的思路和方法与显式R-K法类似。通常，同级的隐式公式获得比显式公式更高的阶。,2020/7/6,68,通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式 R-K法有:,1级2阶中点公式 :,2级2阶梯形公式：,2级4阶R-K公式：,2020/7/6,69,6.变步长方法,在单步法中每一积分步步长实际上是相互独立的，

25、步长的 选择具有了灵活性。根据合理地选择每一积分步的步长， 既保证精度的要求，又可以减少计算量，从而减小舍入误 差。其方便的控制手段是基于误差的事后估计式。,对于给定的精度 ,如果 ,反复将步长折半进行计算，直至 为止，这时再将步长折半一次，就得到所要的结果。,这种通过加倍或折半处理步长的计算方法称为变步长方法。,注：推荐使用精度好计算量低的变步长方法。 用四阶显式R-K方法做变步长方法是实践中较好的方法！,2020/7/6,70,例 分别用改进的欧拉格式和四阶龙格库塔格式解初值问题（取步长h=0.2）：,2020/7/6,71,节点 改进欧拉法 四阶龙格库塔法 准确解 0 1 1 1 0.2

26、 1.186667 1.183229 1.183216 0.4 1.348312 1.341667 1.341641 0.6 1.493704 1.483281 1.483240 0.8 1.627861 1.612514 1.612452 1 1.754205 1.732142 1.732051,表,（注：已指出过准确解,）,2020/7/6,72,单步法的相容性、收敛性和稳定性,1.单步法的相容性 2.单步法的收敛性 3.单步法的稳定性 4.相容性和收敛性的关系 5.相容性和方法阶的关系 6.稳定性和收敛性 7.绝对稳定性和绝对稳定域,2020/7/6,73,单步法的相容性,定义一：对于（

27、9.1.1）常微分方程初值问题 单步法的形式可以变表示为 （9.2.19） 其中 h为步长 若对求解区间中任一固定的x，当 时皆成立,则称由（9.2.19）确定的单步法与微分方程初值问题是相容的 注意到上式左边极限为由（9.1.1）知它应等于从而由相容性 定义得 称相容性条件。,2020/7/6,74,单步法的收敛性,定义二： 设 y(x) 是（9.1.1）的解， 是单步法（9.2.19）产生的数值解，对于每一个固定的 , ， 当 即 。若成立 ， 则称该方法是收敛的。,2020/7/6,75,单步法的稳定性,定义三： 若一个数值方法在基点 处的值有 的扰动，在 此后各基点 （mn）处的值产生

28、的偏差均不超 过 ，则称该方法是稳定的。 单步法的稳定性有以下定理,定理二： 若单步法的增量函数对变量y满足 Lipschtiz 条件， 则单步方法是稳定的。,2020/7/6,76,相容性和收敛性的关系,定理一： 若单步法的增量函数对变量y满足Lipschtiz 条件， 即存在与 h , x 无关的常数 L， 对区域D= 任意两点 （x，y1）,(x,y2)成立，则单步法收敛的充分必要 条件是相容性条件成立。（读者自证）,2020/7/6,77,相容性和方法阶的关系,若单步法是p阶方法则成立 若单步法满足相容性条件，得 所以 =0 也就是说单步法的阶数一定要是正数。由于我们考虑 的单步法皆为

29、正整数，p至少为一。因此我们考虑的 单步法都满足相容性条件。,2020/7/6,78,稳定性和收敛性的关系,若单步法的增量函数满足定理二的条件即单 步法是稳定的则单步法收敛的充分必要条件 是 相容性条件成立。,2020/7/6,79,绝对稳定性和绝对稳定域,稳定性问题是一个比较复杂的问题。为了简化讨论一 般仅对试验方程 进行考察。这里假定 Re0和的允许范围，称为该方法的 绝对稳定域。,2020/7/6,80,绝对稳定性和绝对稳定域2,将Euler方法应用到试验方程得 误差方程是 要求误差不增长则,2020/7/6,81,则Euler 方法的绝对稳定域是以 为半径的圆的内部。为了保证稳定性步长

30、有所控制。 假如当 时h应满足 ，当 时， h 应满 足 等等。 同样我们可以将试验方程用到其它各种单步法当中，求出其绝对稳定域。在实际应用中必须在这个范围内，否则误差传播相当大，即不稳定。,2020/7/6,82,不论单步法或多步法，隐式公式比显式公式稳定性好，但在实际使用隐式公式时，都会遇到两个问题：一个是隐式公式如何能方便地进行计算；另一个是实际计算步长取多大。如隐式梯形公式，每往前推进一步，不必进行多次迭代，而是采用一阶显式Euler公式预测，二阶隐式梯形公式校正一次，构成显式改进Euler公式，能达到与梯形公式同阶的精度，即二阶精度。,3. 预测校正方法,2020/7/6,83,1、数值例题,我们已经学习了很多数值算法，他们的效果到底如何呢？ 下面我们来分析一道例题，看看