【教学课件】第九章常微分方程初值问题数值解法.ppt

《【教学课件】第九章常微分方程初值问题数值解法.ppt》由会员分享，可在线阅读，更多相关《【教学课件】第九章常微分方程初值问题数值解法.ppt（78页珍藏版）》请在三一办公上搜索。

1、第九章 常微分方程初值问题数值解法,9.1 引言9.2 简单的数值方法与基本概念9.3 龙格库塔方法9.4 单步法的收敛性与稳定性9.5 线性多步法9.6 方程组和高阶方程,9.1 引言,本章讨论一阶常微分方程的初值问题：只要函数 适当光滑如满足利普希茨条件：理论上就能保证初值问题的解 存在并且唯一。所谓数值解法，就是寻求解 在一系列离散点上的近似值，相邻两个点间的距离 称为步长。一般情况下我们取为常数，这是节点为：,9.1 引言,初值问题的求解有一个基本特点，它们都是采取“步进式”求解的，即，一步一步地求函数的值。求解的主要方法是：先对方程进行离散化，建立求数值解的递推公式，一类在计算 时只

2、用到前面一步的，称为单步法。另一类在计算 时除了用到理论上就能保证初值问题的解 存在并且唯一。所谓数值解法，就是寻求解 在一系列离散点上的近似值，相邻两个点间的距离 称为步长。一般情况下我们取为常数，这是节点为：,初值问题的求解有一个基本特点，它们都是采取“步进式”求解的，即，一步一步地求函数的值。求解的主要方法是：先对方程进行离散化，建立求数值解的递推公式，一类在计算 时只用到前面一步的，称为单步法。另一类在计算 时除了用到前利用前面一步的，还要用到前面的 等，这种方称为多步法。其次，要研究公式的局部截断误差和阶。数值解和精确解的误差估计和收敛性，还有递推迭代公式的数值稳定性问题。,显式欧拉

3、法 在 平面上，微分方程初值问题的解 称作方程的积分曲线。积分曲线上每一点 的斜率等于该点的函数值。如果按函数 在 平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向一致。方称为多步法。其次，要研究公式的局部截断误差和阶。数值解和精确解的误差估计和收敛性，还有递推迭代公式的数值稳定性问题。,9.2 简单的数值方法与基本概念,显式欧拉法 在 平面上，微分方程初值问题的解 称作方程的积分曲线。积分曲线上每一点 的斜率等于该点的函数值。如果按函数 在 平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向一致。基于上面的几何解释，我们从初始点 出发，先

4、依方向场在该点的方向推进到 上一点，然后再从 依方向场的方向推进到 上一点，循此前进做出一条折线。,9.2 简单的数值方法与基本概念,一般地，设已做出该折线上的顶点，过 依方向场的方向再推进到，显然两个顶点 的坐标有如下关系即 基于上面的几何解释，我们从初始点 出发，先依方向场在该点的方向推进到 上一点，然后再从 依方向场的方向推进到 上一点，循此前进做出一条折线。,一般地，设已做出该折线上的顶点，过 依方向场的方向再推进到，显然两个顶点 的坐标有如下关系即这就是著名的欧拉公式。若初值 已知，则依该公式可逐步算出：,例1 求解初值问题：(其解为)解：根据欧拉方法，得到：这就是著名的欧拉公式。若

5、初值 已知，则依该公式可逐步算出：,例1 求解初值问题：(其解为)解：根据欧拉方法，得到：取步长，计算得到如下结果：,隐式欧拉法 在前面的讨论中，近似计算公式可以看成是由 在区间 上积分得到，而右边的积分是利用左矩形公式 近似，再以 代替 得到，现在右端的积分用右矩形公式，则得到：取步长，计算得到如下结果：,隐式欧拉法 在前面的讨论中，近似计算公式可以看成是由 在区间 上积分得到，而右边的积分是利用左矩形公式 近似，再以 代替 得到，现在右端的积分用右矩形公式，则得到：此式的右端包含，是一种隐式的单步法，称为隐式欧拉方法。利用此法，每一步都要把该式作为 的方程来求解。从数值积分的误差来分析，很

6、难期望隐式欧拉方法比显式欧拉方法精确。为了得到更精确的方法，我们用梯形公式来近似前面的积分，得到梯形方法：它也是一种隐式单步法。,改进的欧拉法 从梯形法的推导，可望它比欧拉法更精确，但它计算量较大，在实际计算中，可取欧拉法的结果为迭代计算 的初值，然后再用梯形公式计算一次，得到：欧拉方法。利用此法，每一步都要把该式作为 的方程来求解。从数值积分的误差来分析，很难期望隐式欧拉方法比显式欧拉方法精确。为了得到更精确的方法，我们用梯形公式来近似前面的积分，得到梯形方法：它也是一种隐式单步法。,改进的欧拉法 从梯形法的推导，可望它比欧拉法更精确，但它计算量较大，在实际计算中，可取欧拉法的结果为迭代计算

7、 的初值，然后再用梯形公式计算一次，得到：或写成：改进的欧拉法是一种显式单步法，有时为了计算方便，也用下面的公式：其中：,单步法的截断误差与阶 初值问题的单步法可用一般形式表示为：其中多元函数 与 有关，当 含有 时，方法是隐式的，不含 时则为显式的，所以显式单步法可表示为：称 为增量函数，对于欧拉法，或写成：改进的欧拉法是一种显式单步法，有时为了计算方便，也用下面的公式：其中：,单步法的截断误差与阶 初值问题的单步法可用一般形式表示为：其中多元函数 与 有关，当 含有 时，方法是隐式的，不含 时则为显式的，所以显式单步法可表示为：称 为增量函数，对于欧拉法，为了分析其截断误差，我们将 在 处

8、作泰勒展开，得到：假设 是精确的，于是可得欧拉法计算公式的误差为：,定义1 设 是初值问题的准确解，称为显式单步法 的局部截断误差。之所以称为局部的，是假设在 前各步没有误差，当 时，计算一步，则有开，得到：假设 是精确的，于是可得欧拉法计算公式的误差为：,定义1 设 是初值问题的准确解，称为显式单步法 的局部截断误差。之所以称为局部的，是假设在 前各步没有误差，当 时，计算一步，则有所以，局部截断误差可以理解为用计算的方法计算一步的误差，这样，欧拉法的局部截断误差为：这里 称为局部截断误差主项，显然 一般情况下的定义如下：,定义2 设 是初值问题的准确解，若存在最大整数 使显式单步法的局部误

9、差满足则称该方法具有 阶精度。若将上式展开成则 称为局部截断误差主项。所以，局部截断误差可以理解为用计算的方法计算一步的误差，这样，欧拉法的局部截断误差为：这里 称为局部截断误差主项，显然 一般情况下的定义如下：,定义2 设 是初值问题的准确解，若存在最大整数 使显式单步法的局部误差满足则称该方法具有 阶精度。若将上式展开成则 称为局部截断误差主项。以上定义对隐式单步法也适用，如对后退欧拉法，其局部截断误差为：这里，是1阶方法，局部截断误差主项为,同样，对梯形方法有：所以梯形方法是2阶的，其局部误差主项为：以上定义对隐式单步法也适用，如对后退欧拉法，其局部截断误差为：这里，是1阶方法，局部截断

10、误差主项为,同样，对梯形方法有：所以梯形方法是2阶的，其局部误差主项为：改进欧拉法的例 梯形法精度得到提高，但算法复杂，计算量大。为了控制计算量，通常只迭代一两次就转入下一步的计算，这样就可以简化计算。具体做法：先用欧拉公式求一个初值，称为预测值，再用梯形公式校正一次，这样建立的公式通常称为改进的欧拉公式：,预测：校正：或表示为：所以梯形方法是2阶的，其局部误差主项为：改进欧拉法的例 梯形法精度得到提高，但算法复杂，计算量大。为了控制计算量，通常只迭代一两次就转入下一步的计算，这样就可以简化计算。具体做法：先用欧拉公式求一个初值，称为预测值，再用梯形公式校正一次，这样建立的公式通常称为改进的欧

11、拉公式：,预测：校正：或表示为：例：用改进的欧拉方法求解初值问题：解：仍取，计算结果如表,显式龙格库塔法的一般形式 欧拉法的局部截断误差为，是1阶的方法对于改进的欧拉法，它可表示为：此时增量函数例：用改进的欧拉方法求解初值问题：解：仍取，计算结果如表,9.3 龙格库塔方法,显式龙格库塔法的一般形式 欧拉法的局部截断误差为，是1阶的方法对于改进的欧拉法，它可表示为：此时增量函数它比欧拉法的 多计算了一个函数值，可望，若要使得到的公式阶数更大，就必须包含更多的 值，实际上从方程 等价的积分形式，即(*),9.3 龙格库塔方法,要使公式的阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加

12、求积节点，为此可将(*)的右端用求积公式表示为：一般来说，点数 越多，精度越高，上式右端相当于增量函数，为得到便于计算的显式方法，可它比欧拉法的 多计算了一个函数值，可望，若要使得到的公式阶数更大，就必须包含更多的 值，实际上从方程 等价的积分形式，即(*),要使公式的阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积节点，为此可将(*)的右端用求积公式表示为：一般来说，点数 越多，精度越高，上式右端相当于增量函数，为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为：其中这里 均为常数，方法称为 级显式R-K方法,当 时，就是欧拉法，此时 当 时，改进欧拉法就是其中的一

13、种，下面证明，要使公式有更高的阶，就要增加点数。下面我们仅就 推导R-K方法，并给出 时的常用公式，其推导方法与 时类似，只是计算较为复杂。增量函数，为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为：其中这里 均为常数，方法称为 级显式R-K方法,当 时，就是欧拉法，此时 当 时，改进欧拉法就是其中的一种，下面证明，要使公式有更高的阶，就要增加点数。下面我们仅就 推导R-K方法，并给出 时的常用公式，其推导方法与 时类似，只是计算较为复杂。二阶显式R-K方法 对于 的R-K方法，可以得到如下的计算公式：这里 均为待定常数，我们希望适当选取这些系数，使公式阶数 尽量高。根据局部截断误差

14、的定义，上式的局部截断误差为：,这里 为得到 的阶，要将上式各项在 处做泰勒展开，由于 是二元函数，故要用到二元泰勒展开，各项展开为：二阶显式R-K方法 对于 的R-K方法，可以得到如下的计算公式：这里 均为待定常数，我们希望适当选取这些系数，使公式阶数 尽量高。根据局部截断误差的定义，上式的局部截断误差为：,这里 为得到 的阶，要将上式各项在 处做泰勒展开，由于 是二元函数，故要用到二元泰勒展开，各项展开为：其中,将以上结果代入：得到：其中,将以上结果代入：得到：要使公式成为2阶的，必须有：,即：上式的解不唯一，可以令，则得这样得到的公式称为二阶R-K方法，如取，则，这就是改进的欧拉法。若取

15、，则，得到：要使公式成为2阶的，必须有：,即：上式的解不唯一，可以令，则得这样得到的公式称为二阶R-K方法，如取，则，这就是改进的欧拉法。若取，则，得到：称为中点公式，相当于数值积分的中矩形公式，也即,对 的R-K公式能否使局部误差提高到？为此需要把 多展开一项，通过分析可以得到不可能的结论。因此 时只能得到2阶的公式。三阶与四阶显式R-K方法 要得到三阶的显式R-K方法，必须，此时计算 若取，则，得到：称为中点公式，相当于数值积分的中矩形公式，也即,对 的R-K公式能否使局部误差提高到？为此需要把 多展开一项，通过分析可以得到不可能的结论。因此 时只能得到2阶的公式。三阶与四阶显式R-K方法

16、 要得到三阶的显式R-K方法，必须，此时计算公式为：其中 及 均为待定参数，误差为,只要将 按二元函数的泰勒公式展开，使 可得待定参数需满足的方程为：要得到三阶的显式R-K方法，必须，此时计算公式为：其中 及 均为待定参数，误差为,只要将 按二元函数的泰勒公式展开，使 可得待定参数需满足的方程为：这里为8个未知数，6个方程，解也是不唯一的，可以得到很多公式，下面是其中称为库塔的公式,这里为8个未知数，6个方程，解也是不唯一的，可以得到很多公式，下面是其中称为库塔的公式,继续上述过程，经过较复杂的演算，可以导出各种四阶龙格库塔公式，下面是常用的经典公式：,例3：用四阶龙格库塔方法求解定解问题：取

17、步长，从 到。解：写出具体的四阶龙格库塔的具体格式为：继续上述过程，经过较复杂的演算，可以导出各种四阶龙格库塔公式，下面是常用的经典公式：,例3：用四阶龙格库塔方法求解定解问题：取步长，从 到。解：写出具体的四阶龙格库塔的具体格式为：,列出计算结果如下表：,列出计算结果如下表：可见计算精度较以前大为提高。变步长的龙格库塔方法 单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在求解区间内的步数就增加了，这在增加计算量的同时也增加了计算的舍入误差，因此也有一个选择步长的问题。,选择步长需要考虑的两个问题：怎样衡量和检验计算结果的精度？如何依据所获得的精度处理步长？我们考虑经典的四阶龙格库

18、塔公式，从节点 出发，先以 为步长求出一个近似值，记为，由于公式的局部截断误差为，故变步长的龙格库塔方法 单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在求解区间内的步数就增加了，这在增加计算量的同时也增加了计算的舍入误差，因此也有一个选择步长的问题。,选择步长需要考虑的两个问题：怎样衡量和检验计算结果的精度？如何依据所获得的精度处理步长？我们考虑经典的四阶龙格库塔公式，从节点 出发，先以 为步长求出一个近似值，记为，由于公式的局部截断误差为，故然后将步长折半，即取 为步长，从 跨两步到再来求得一个近似值，每一步的截断误差是因此有,比较前后两式，我们得到，步长折半后，误差大约减少到

19、，即有由此易得下列事后估计式：这样，我们可以通过检查步长，折半前后两次计算结果的偏差 来判定所选的步长是否合适，公式的局部截断误差为，故然后将步长折半，即取 为步长，从 跨两步到再来求得一个近似值，每一步的截断误差是因此有,比较前后两式，我们得到，步长折半后，误差大约减少到，即有由此易得下列事后估计式：这样，我们可以通过检查步长，折半前后两次计算结果的偏差 来判定所选的步长是否合适，分下列两种情况处理，称为变步长方法：1)对于给定的精度，如果，我们反复将步长折半进行计算，直到 为止，这时取最终得到的作为结果；2)如果，我们将反复将步长加倍，直到 为止，这时再将步长折半一次，得到结果。,9.4

20、单步法的收敛性与稳定性,收敛性与相容性 数值解法的基本思想是通过离散将微分方程转化为差分方程，如单步法(1)它在 处的解为，而初值问题的解为，记分下列两种情况处理，称为变步长方法：1)对于给定的精度，如果，我们反复将步长折半进行计算，直到 为止，这时取最终得到的作为结果；2)如果，我们将反复将步长加倍，直到 为止，这时再将步长折半一次，得到结果。,9.4 单步法的收敛性与稳定性,收敛性与相容性 数值解法的基本思想是通过离散将微分方程转化为差分方程，如单步法(1)它在 处的解为，而初值问题的解为，记，称为整体截断误差，收敛性就是讨论当 固定，且 时 的问题。定义3 若一种数值方法对于固定的，当时

21、有，则称该方法是收敛的。显然，数值方法收敛是指，对单步法(1)有下述收敛性定理：,9.4 单步法的收敛性与稳定性,定理1 假设单步法(1)具有 阶精度，且增量函数 关于 满足利普希茨条件(2)又设初值 准确，即，则其整体截断误差为(3)证明：设以 表示取 用公式(1)求得的结果,固定，且 时 的问题。定义3 若一种数值方法对于固定的，当时有，则称该方法是收敛的。显然，数值方法收敛是指，对单步法(1)有下述收敛性定理：,定理1 假设单步法(1)具有 阶精度，且增量函数 关于 满足利普希茨条件(2)又设初值 准确，即，则其整体截断误差为(3)证明：设以 表示取 用公式(1)求得的结果,即(4)则

22、为局部截断误差，由于所给的方法具有 阶精度，按定义2，存在定数，使又由(4)和(1)，有,利用条件(2)，有：从而有即对整体截断误差 成立下列递推关系：反复递推，可得：证明：设以 表示取 用公式(1)求得的结果,即(4)则 为局部截断误差，由于所给的方法具有 阶精度，按定义2，存在定数，使又由(4)和(1)，有,利用条件(2)，有：从而有即对整体截断误差 成立下列递推关系：反复递推，可得：再注意到当 时，最终得到下列估计式：由此可以断定，如果初值是准确的，则(3)成立。,推论：对于一个 阶的显式单步法，若微分方程的右端函数 关于 满足利普希茨条件，且初值是精确的，则显式欧拉法，改进欧拉法和龙格

23、库塔方法是收敛的。定理1表明，时单步法收敛，并且当 是初值问题的解，且具有 阶精度时，有展开式：再注意到当 时，最终得到下列估计式：由此可以断定，如果初值是准确的，则(3)成立。,推论：对于一个 阶的显式单步法，若微分方程的右端函数 关于 满足利普希茨条件，且初值是精确的，则显式欧拉法，改进欧拉法和龙格库塔方法是收敛的。定理1表明，时单步法收敛，并且当 是初值问题的解，且具有 阶精度时，有展开式：所以 的充要条件是，而，于是可给出如下的定义：定义4 单步法的增量函数 称为该单步法与初值问题相容。,以上讨论表明 阶方法当 与初值问题相容，反之，相容的方法至少是一阶的。于是，由定理1可知，线性单步

24、方法收敛的充分必要条件是该方法是相容的。绝对稳定性与绝对稳定域 所以 的充要条件是，而，于是可给出如下的定义：定义4 单步法的增量函数 称为该单步法与初值问题相容。,以上讨论表明 阶方法当 与初值问题相容，反之，相容的方法至少是一阶的。于是，由定理1可知，线性单步方法收敛的充分必要条件是该方法是相容的。绝对稳定性与绝对稳定域 定义5 若一种数值方法在节点值 上有大小为 的扰动，而于以后各节点上产生的偏差均不超过，则称该方法是稳定的。例4 考察初值问题,其准确解 是一个按指数曲线衰减的很快的函数，如图：绝对稳定性与绝对稳定域 定义5 若一种数值方法在节点值 上有大小为 的扰动，而于以后各节点上产

25、生的偏差均不超过，则称该方法是稳定的。例4 考察初值问题,其准确解 是一个按指数曲线衰减的很快的函数，如图：,用欧拉方法解该方程，得到：若取，则欧拉法的具体形式为：若取，则欧拉法的具体形式为：显然，前一形式计算不稳定，后一形式计算稳定。再考察后退欧拉方法，可以得到，当 时，计算是稳定的。具体数据见书上356页表94。例题表明：稳定性不但与方法有关，也与步长 有关。当然也与 有关，为了只考察数值方法本身，通常只检验将数值方法用于解模型方程的稳定性，模型方程为：，其中 为复数。,对于一般的方程，将方程 中的 在解域内某点 作泰勒展开并局部线性化：忽略高阶项，令 和 对上式显然，前一形式计算不稳定，

26、后一形式计算稳定。再考察后退欧拉方法，可以得到，当 时，计算是稳定的。具体数据见书上356页表94。例题表明：稳定性不但与方法有关，也与步长 有关。当然也与 有关，为了只考察数值方法本身，通常只检验将数值方法用于解模型方程的稳定性，模型方程为：，其中 为复数。,对于一般的方程，将方程 中的 在解域内某点 作泰勒展开并局部线性化：忽略高阶项，令 和 对上式作变量代换，得到：这就是模型方程。首先研究欧拉方程的稳定性。模型方程的欧拉公式为：假设在节点 上有一扰动，它的传播使节点值 产,生大小为 的扰动值，假设用 按欧拉公式得出 的计算过程不再有新的误差，则扰动值满足可见扰动值满足原来的差分方程，这样

27、如果差分方程的解是不增长的，即有：这就是模型方程。首先研究欧拉方程的稳定性。模型方程的欧拉公式为：假设在节点 上有一扰动，它的传播使节点值 产,生大小为 的扰动值，假设用 按欧拉公式得出 的计算过程不再有新的误差，则扰动值满足可见扰动值满足原来的差分方程，这样如果差分方程的解是不增长的，即有：则它就是稳定的。这一结论对于以下研究的其它方法也适用。显然，为了保证，只要，这在的复平面上是一个以(-1,0)为圆心，1为半径的单位圆我们称其为欧拉法的绝对稳定域，一般情形可有如下的定义。,定义6 单步法 用于模型，若得到的解，满足，则称该方法是绝对稳定的。使得 成立的区域，称为绝对稳定区域，它与实轴的交

28、称为绝对稳定区间。对于欧拉法，其绝对稳定域已由给出，绝对稳定区间为。在具体问题中，由则它就是稳定的。这一结论对于以下研究的其它方法也适用。显然，为了保证，只要，这在的复平面上是一个以(-1,0)为圆心，1为半径的单位圆我们称其为欧拉法的绝对稳定域，一般情形可有如下的定义。,定义6 单步法 用于模型，若得到的解，满足，则称该方法是绝对稳定的。使得 成立的区域，称为绝对稳定区域，它与实轴的交称为绝对稳定区间。对于欧拉法，其绝对稳定域已由给出，绝对稳定区间为。在具体问题中，由于 不同，因此，是否稳定就与 有关。对于二阶龙格库塔方法，解模型方程可得：故：,绝对稳定域由 得到，于是可得到为即，类似可得三

29、阶及四阶的龙格库塔方法的 分别为：于 不同，因此，是否稳定就与 有关。对于二阶龙格库塔方法，解模型方程可得：故：,绝对稳定域由 得到，于是可得到为即，类似可得三阶及四阶的龙格库塔方法的 分别为：由 可得相应的绝对稳定域，从而得到绝对稳定区间如下：三阶龙格库塔方法：四阶龙格库塔方法：龙格库塔方法的稳定域为有限域，都对步长有限制,自学书上359页的例，学会稳定区间的具体求法。9.5 线性多步法 不作要求。由 可得相应的绝对稳定域，从而得到绝对稳定区间如下：三阶龙格库塔方法：四阶龙格库塔方法：龙格库塔方法的稳定域为有限域，都对步长有限制,自学书上359页的例，学会稳定区间的具体求法。9.5 线性多步

30、法 不作要求。,9.6 方程组和高阶方程,一阶方程组 前面关于一阶方程的所有解法，都可用于一阶方程组，这只要把 和 理解为向量就可以了。考察一阶方程组：初始条件为：这时，龙格库塔法为：其中：,9.6 方程组和高阶方程,自学书上375页两个方程的情形，加深理解方程组的解法。高阶方程 化为一阶方程组来求解考虑高阶方程：初始条件为：这时，龙格库塔法为：其中：,自学书上375页两个方程的情形，加深理解方程组的解法。高阶方程 化为一阶方程组来求解考虑高阶方程：初始条件为：只要引进新变量就可以把原来的高阶方程化为一阶方程组,此时，初始条件也相应的化为：根据微分方程的理论，我们知道它们是等价的，因此问题得到解决。只要引进新变量就可以把原来的高阶方程化为一阶方程组,此时，初始条件也相应的化为：根据微分方程的理论，我们知道它们是等价的，因此问题得到解决。自学书上376页关于二阶方程的叙述，理解高阶方程的求解过程。刚性方程部分不作要求。,