常微分方程初值问题的数值解法

常微分方程初值问题的数值解法

第6章1可编辑ppt引言在实际问题中，常需要求解微分方程(如发电机转子运动方程)。只有简单的和典型的微分方程可以求出解析解，而在实际问题中的微分方程往往无法求出解析解。常微分方程：-----------(1)-----------(2)一阶常微分方程2可编辑ppt-----------(3)(1),(2)式称为初值问题,(3)式称为边值问题-----------(4)另外,在实际应用中还经常需要求解常微分方程组:

本课程主要研究问题一阶常微分方程(1)的数值解法,我们首先介绍初值问题(1)的解存在的条件3可编辑ppt定理只要f(x,y)连续，且关于y

满足Lipschitz

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

使对任意定义在[a,b]上的y1(x)和y2(x)都成立，则初值问题（1）存在唯一解。（通常采用等距节点）对于问题(1)要求它的数值解4可编辑ppt常微分方程数值解公式的推导

求初值问题数值解的方法是步进法，即从已知的初值y0出发，通过一定的计算求y1

，然后由y1或y0和y1求出y2

，依次计算到yn

，即在计算出yk后计算yk+1

，这时有单步法：计算yk+1时，只利用yk多步法：计算yk+1时，用到yk,yk-1,yk-2,…常微分方程数值解公式的主要推导方法泰勒展开利用差商利用数值积分法5可编辑ppt1、泰勒展开的求解思路：将按泰勒级数展开用的近似值代入上式右端，记所得结果为，则得到数值解序列的计算公式：6可编辑ppt2、化导数为差商的求解方法思路：若在点处的导数用差商来近似代替，如向前差商则微分方程初值问题化为将近似号改为等号，精确解

改为近似解，得7可编辑ppt3、数值积分的求解思路：如果将微分方程在各小区间上对其两边进行积分，即如用矩形数值积分公式可得：8可编辑ppt以上三种方法推导出同一个数值求解公式:这个数值公式称为欧拉(Euler)公式。9可编辑ppt§6.1欧拉方法一、

欧拉格式：x0x1向前差商近似导数记为欧拉公式几何意义

用一条通过初始点的折线近似表示解曲线

,亦称为欧拉折线法，或称为矩形法。一般形式1、显式欧拉公式10可编辑ppt11可编辑ppt在假设yk=y(xk)，即第

k

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

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

阶精度。定义

欧拉法的局部截断误差：欧拉法具有1阶精度。局部截断误差和阶数12可编辑ppt2、隐式欧拉格式向后差商近似导数x0x1))(,()(1101xyxfhyxy+)1,...,0(),(111-=+=+++nkyxfhyykiik由于未知数yk+1

同时出现在等式的两边，不能直接得到，故称为隐式欧拉公式，而前者称为显式欧拉公式。一般先用显式计算一个初值，再迭代求解。

隐式欧拉法的局部截断误差：即隐式欧拉公式具有1阶精度。13可编辑ppt二、两步欧拉格式（中点公式）中心差商近似导数x0x2x1假设,则可以导出即两步欧拉格式具有2阶精度。该方法需要2个初值y0和y1来启动递推过程，这样的算法称为双步法。14可编辑ppt三、

梯形公式—显、隐式两种算法的平均注：有局部截断误差，即梯形公式具有2

阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，不易求解。

对欧拉法进行改进，用梯形公式计算右侧积分，即计算公式15可编辑ppt梯形格式算法计算步骤：先用(1)式计算出处。再用(2)式反复进行迭代，得到计算公式-----------(1)-----------(2)类似地得到用控制迭代次数，为允许误差。把满足误差要求的作为的近似值。16可编辑ppt方法显式欧拉法隐式欧拉法梯形公式中点公式简单精度低稳定性最好精度低,计算量大精度提高计算量大精度提高,显式多一个初值,可能影响精度不同方法比较17可编辑ppt四、改进欧拉法（预报-校正法）Step1:

先用显式欧拉公式作预报，算出),(1iiiiyxfhyy+=+Step2:再将代入隐式梯形公式的右边作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy它可表示为嵌套形式表示为平均化形式此法称为预报-校正法，是显式算法。注：可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式（只迭代一次）

，比隐式梯形公式的迭代求解过程简单。脚标用i18可编辑ppt举例：19可编辑ppt20可编辑pptxEuler法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.4832400.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.7378671.73205121可编辑ppt§6.2龙格-库塔法一、泰勒级数法

龙格—库塔(Runge-Kutta)法(简称为R-K方法)是一类高精度的一步法，这类方法与泰勒级数法有着密切的关系。

设有初值问题由泰勒展开式

从理论上讲，只要解y(x)有任意阶导数，泰勒展开方法就可以构造任意阶求yk+1公式。但由于计算这些导数是非常复杂的，所以这种方法实际上不能用来解初值问题。

22可编辑ppt设有初值问题二、龙格－库塔法的基本思路等价于：（积分中值定理）

R-K方法基本思想：用在几个不同点的加权平均值（线性组合）来代替准确的的值，构造近似公式。再把近似公式与解的泰勒展开式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。这样龙格－库塔法保留了泰勒级数展开法的高阶局部截断误差，又避免了高阶导数的计算。我们先分析欧拉法与预估—校正法。23可编辑ppt24可编辑ppt推广这种单步法称为Runge-Kutta方法,简记为R-K公式.Ki为某些点上的斜率，或f(x,y)在某些点上的值。25可编辑ppt三、二阶龙格-库塔法目标：建立高精度的单步递推格式。单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1

,yi+1

)点。欧拉法及其各种变形所能达到的最高精度为2阶。

考察改进的欧拉法，可以将其改写为：斜率一定取K1K2的平均值吗？步长一定是一个h

吗？脚标用i26可编辑ppt首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得

Step1:将K2在(xi,yi)

点作Taylor展开将改进欧拉法推广为：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii++==++=+llStep2:将K2代入yi+1表达式，得到27可编辑pptStep3:将yi+1与y(xi+1)在xi点的泰勒展开作比较要求，则必须有：这里有个未知数，个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。注意到，就是改进的欧拉法。Q:

为获得更高的精度，应该如何进一步推广？28可编辑ppt其中i

(i=1,…,m)，i

(i=2,…,m)

和ij

(i=2,…,m;j=1,…,i1

)

均为待定系数，确定这些系数的步骤与前面相似。其解不唯一。)...,(......),(),(),(]...[1122112321313312122122111--++++++=+++=++==++++=mmmmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyybbbabbaballl最常用为四阶经典龙格-库塔法29可编辑ppt四阶经典龙格-库塔法公式四、四阶龙格-库塔法用四个f函数值的线性组合得到四阶龙格-库塔法。经典龙格-库塔法公式具有四阶精度，因此可取大步长。30可编辑ppt注：

龙格-库塔法的主要运算在于计算Ki

的值，即计算f

的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：753可达到的最高精度642每步须算Ki的个数高于四阶时每步计算量增加较多，但精度提高不快，因此使用的比较少。由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h

取小。31可编辑ppt§3收敛性与稳定性/*ConvergencyandStability*/

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

(同时i)时有yi

y(xi

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

，有32可编辑ppt稳定性/*Stability*/例：考察初值问题在区间[0,0.5]上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解改进欧拉法

欧拉隐式欧拉显式

节点xi

1.00002.00004.00008.00001.6000101

3.2000101

1.00002.500010