常微分方程的数值解法

常微分方程的数值解法1第一页，共三十三页，2022年，8月28日上一页下一页返回

本章介绍求解微分方程数值解的基本思想和方法.

含有自变量、未知函数和它的一阶导数和高阶导数的方程.常微分方程它是描述运动、变化规律的重要数学方法之一，分为两类：1.初值问题即给出未知函数及导数在初始点的值；2.边值问题即给出未知函数及（或）它的某些导数在区间两个端点的值。2第二页，共三十三页，2022年，8月28日

考虑一阶常微分方程的初值问题:只要f(x,y)在[a,b]R1上连续，且关于y

满足Lipschitz

条件，即存在与x,y无关的常数L

使对任意定义在[a,b]上的y1(x)和y2(x)都成立，则上述问题解存在唯一解。所谓数值解法就是要计算出初值问题的解函数y(x)在一系列离散点a=x0<x1<…<xN=b上的近似值：y0,y1,……yN.节点间距为步长，通常采用等距节点，即取hi=h

(常数)。{yn}称为问题的数值解.数值解所满足的离散方程统称为差分格式.

上一页下一页返回

3第三页，共三十三页，2022年，8月28日第一节

欧拉方法一、欧拉公式令yn为y(xn)的近似值，将上式代入(*)式可得此式称为欧拉(Euler)公式.

为Euler方法的局部截断误差.上一页下一页返回

4第四页，共三十三页，2022年，8月28日例1

用欧拉公式解初值问题解：取步长h=0.1,欧拉公式的具体形式为：依次计算可得………

y3，y4，y5，y6，y7，y8，y9，y10上一页下一页返回

5第五页，共三十三页，2022年，8月28日其部分结果见下表

可见Euler方法的计算结果精度不太高。

上一页下一页返回

6第六页，共三十三页，2022年，8月28日欧拉公式的几何意义：x0P0x1P1x2P2xnPn几何意义：用折线近似代替方程的解曲线,因而也称Euler方法为折线法.上一页下一页返回

7第七页，共三十三页，2022年，8月28日二、后退的欧拉公式也用一阶差商逼近导数令yn+1为y(xn+1)的近似值，则可得称为后退Euler公式已知yn时,必须通过解方程才能求出yn＋1

,这样的公式称为隐式公式,而Euler公式为显式公式.

Euler公式和后退Euler公式都是由yn去计算yn+1，因此，称它们为单步法。上一页下一页返回

8第八页，共三十三页，2022年，8月28日定义在假设yi=y(xi)，即第

i

步计算是精确的前提下，考虑的截断误差Ti+1=y(xi+1)

yi+1称为局部截断误差。定义若某算法的局部截断误差为O(hp+1)，则称该算法有p

阶精度。显然,p越大,精度越高.三、局部截断误差与方法的阶（将准确解代入公式的左、右两端，其左端与右端之差）

Euler方法的精度

其中：

上一页下一页返回

9第九页，共三十三页，2022年，8月28日所以，Euler方法具有1阶精度。将在点处一阶Taylor展开上一页下一页返回

10第十页，共三十三页，2022年，8月28日所以，后退的Euler方法也具有1阶精度。将在点处一阶Taylor展开隐式Euler方法的精度

上一页下一页返回

11第十一页，共三十三页，2022年，8月28日—显、隐式两种算法的平均

欧拉公式的改进其局部误差为：此公式具有2阶精度.称平均公式或梯形公式梯形公式可由下迭代式计算：其中迭代初值是Euler公式提供.上一页下一页返回

12第十二页，共三十三页，2022年，8月28日四、改进的欧拉公式Step1:

先用显式欧拉公式作预测，算出),(1iiiiyxfhyy+=+Step2:再将代入隐式梯形公式的右边作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy注：此法亦称为预测-校正法。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。它的精度高于显式欧拉法。上一页下一页返回

13第十三页，共三十三页，2022年，8月28日为了便于编程,常将改进的欧拉公式写为：上一页下一页返回

14第十四页，共三十三页，2022年，8月28日例2用改进的欧拉法解例1中的初值问题.解：取步长h=0.1,

改进欧拉法的具体形式为具体计算过程如下上一页下一页返回

15第十五页，共三十三页，2022年，8月28日xn改进的欧拉法误差xn改进的欧拉法误差0100.61.4859560.0027160.21.1840960.0000880.81.6164760.0040240.41.3433600.0017191.01.7378690.005818依次计算可得………

y3，y4，y5，y6，y7，y8，y9，y10其部分结果见下表

上一页下一页返回

16第十六页，共三十三页，2022年，8月28日例3对下面的初值问题解

（1）取步长h=0.1,欧拉方法的具体公式为（2）取步长h=0.1,改进的欧拉方法的具体公式为取步长h=0.1，分别用Euler方法、改进的Euler方法求数值解。上一页下一页返回

17第十七页，共三十三页，2022年，8月28日计算结果见下表Euler方法改进的Euler方法xnynyn0.10.20.30.40.50.60.70.80.91.00.9000000.8100000.7290000.6561000.5904900.5314410.4782970.4304670.3874210.3486790.9050000.8190250.7412180.6708020.6070760.5494040.4972100.4499750.4072280.368541上一页下一页返回

18第十八页，共三十三页，2022年，8月28日第二节

龙格-库塔法基本思想考察改进的欧拉法，可以将其改写为：斜率一定取k1k2的平均值吗？步长一定是一个h

吗？只要能对平均斜率提供一种近似算法,就能得到一种对应的差分格式.上一页下一页返回

19第十九页，共三十三页，2022年，8月28日例如取m个点的斜率构造如下形式的公式该公式称为m级龙格－库塔(Runge-Kutta)公式,简称R-K公式.求解：只需将公式的局部截断误差在xn点进行Taylor展开,令其前面尽可能多的项为0,便可导出ai，bij，ci所满足的方程组,即可从中求出这些系数.

上一页下一页返回

20第二十页，共三十三页，2022年，8月28日以m=2的情形为例说明建立R-K公式的方法.其局部截断误差为：上一页下一页返回

21第二十一页，共三十三页，2022年，8月28日因此有：而对于h3，若将k2的Taylor展开式多取一项,会发现h3项的系数不可能为0.而对于上式有无穷多个解,它的每一组解都给出了一个局部截断误差为的二级R-K公式,即二阶R-K公式.当取时,二阶R-K公式就是改进的Euler公式

这里有个未知数，个方程。32上一页下一页返回

22第二十二页，共三十三页，2022年，8月28日常用的标准四阶R－K公式（经典R-K方法）最常用的四阶标准R－K公式（经典R-K方法）为：上一页下一页返回

23第二十三页，共三十三页，2022年，8月28日例用四阶标准R-K公式解初值问题

解：取h=0.2，四阶标准R-K法的具体格式如下:上一页下一页返回

24第二十四页，共三十三页，2022年，8月28日已知

上一页下一页返回

25第二十五页，共三十三页，2022年，8月28日同理可计算得具体结果见下表至少具有四位有效数字.比较:上节用改进的Euler公式计算,取h=0.1，最多具有四位有效数字

。

上一页下一页返回

26第二十六页，共三十三页，2022年，8月28日改进的Euler公式每前进一步只要计算两次f值,而4阶R-K公式每前进一步要计算四次f值,但改进的Euler法的步长比4阶R-K法的小一半,两者计算总量差不多.

而4阶R-K法的效果要比改进的Euler法好.

由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h

取小。上一页下一页返回

27第二十七页，共三十三页，2022年，8月28日第三节

单步法的收敛性与稳定性收敛性

/*Convergency*/定义若某算法对于任意固定的x=xi=x0+ih，当h0

(同时i)时有yi

y(xi

)，则称该算法是收敛的。例：就初值问题考察欧拉显式格式的收敛性。解：该问题的精确解为欧拉公式为对任意固定的x=xi=ih

，有上一页下一页返回

28第二十八页，共三十三页，2022年，8月28日稳定性

/*Stability*/例：考察初值问题在区间[0,0.5]上的解。分别用欧拉法、隐式欧拉法和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解改进欧拉法隐式欧拉法欧拉法

节点xi

1.00002.00004.00008.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107Whatiswrong??!上一页下一页返回

29第二十九页，共三十三页，2022年，8月28日定义若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的.一般分析时为简单起见，只考虑试验方程常数l<0，可以是复数当步长取为h

时，将某算法应用于上式，并假设在初值产生误差，则若此误差以后逐步衰减，就称该算法相对于z=lh

绝对稳定，z

的全体构成绝对稳定区域。我们称算法A

比算法B

稳定，就是指A的绝对稳定区域比B