郑州大学研究生课程数值分析复习---第八章 常微分方程数值解法
08常微分方程数值解法(工程数学)
k2为 x i + ph 点处的切线斜率值,对照改进的欧拉 法,将
x i + p 视为 x i +1，即可得 k 2 = f ( x i + ph, y i + phk1 )
21
对常微分方程初值问题的解 y=y(x),根据微分中值定 理，存在点 ξ ∈ ( x i , x i + 1 ) ，使得
0.96 0.8523 0.7003 0.5325 0.3750
0.9616 0.8549 0.7025 0.5327 0.3725
0.9615 0.8549 0.7022 0.5327 0.3730
0.9608 0.8521 0.6977 0.5273 0.3679
3.1、龙格-库塔法
考虑用函数f(x,y)在若干点上的函数值的线性组合 来构造近似公式，构造时要求近似公式在(xi,yi)处的 Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项 重合，从而使近似公式达到所需要的阶数。既避免求偏 导,又提高了计算方法精度的阶数。或者说,在[xi,xi+1]这 一步内多预报几个点的斜率值，然后将其加权平均作为 平均斜率，则可构造出更高精度的计算格式，这就是龙 格—库塔（Runge-Kutta）法的基本思想。
20
1 二阶龙格—库塔法 在 [x , x ] 上取两点xi和 xi+ p = xi + ph,以该两点处
i i +1
的斜率值k1和k2的加权平均(或称为线性组合)来求取 平均斜率k*的近似值K，即
K = λ 1k 1 + λ 2k 2
式中:k1为xi点处的切线斜率值，k = f ( x , y ) = y ′ ( x )
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
常微分方程数值解法
介绍常微分方程数值解法常微分方程（ordinary differential equations，ODE）可用于描述许多日常存在的物理系统。
处理ODE问题常常被称为数值求解法，这指的是找到概括ODE或者其他适用于数学模型的解决方案来模括这些ODE。
这种解决方案可能在一系列不同方案中发挥重要作用，以此来提供更好的解释和预测。
常微分方程与几何图形更为相关，它利用二维或者三维空间中曲线的绘制以及分析。
通过引入一些不同的方法，可以对不同的常微分方程中的量进行描述，使得可以通过数值方法的解析来进行研究。
数值解法可能是时间消耗较多的，但有助于验证几何图形中的某些过程，以此帮助揭示数学模型。
四种常见的常微分方程数值解法四种常见的常微分方程数值解法是：前向差分法、向后差分法、中点法和全分方法。
•前向差分法：前向差分法的基本概念是利用ODE的特定解来表达时间步的影响。
这是一种基本的数值法，可以在ODE中确定任意位置的点作为终点。
在这里，任何这样的点都可以表示为ODE右边的时间步。
•向后差分法：它是反过来基于前向差分法。
它要求对ODE中的时间步进行逆向推导，以获得某一特定点的解。
向后差分法要求推导反向解中点，以便可以从每一步中获取该点的解。
•中点法：这是一种非常基本的数值解法，可以用来求解ODE中的某一步的解，但不具有直观的方法解释。
主要的思想是在每一次时间步中通过求出ODE的中点来寻找解。
•全分方法：这是一种更复杂的数值解法，它要求将ODE中的每一步解细分并解决。
与前面提到的三种解法不同，它首先要求将ODE分解成若干离散区间，然后计算每一段区间中的点。
这种解法可以更准确地进行处理，但时间消耗较多，因此比较少被使用。
优化方案在需要解决常微分方程时，为了得到最佳的结果，有必要考虑一些优化措施。
•首先，应考虑将一个复杂的ODE拆分成一些更易解决的问题。
这样做的结果是，预见到解决此ODR的总耗时将会降低。
•其次，为了加快计算速度，可以考虑使用预解算法。
常微分方程数值解法
i.常微分方程初值问题数值解法i.1 常微分方程差分法考虑常微分方程初值问题：求函数()u t 满足(,), 0du f t u t T dt=<≤ （i.1a ） 0(0)u u = (i.1b)其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。
我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。
通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。
本章讨论常微分方程最常用的近似数值解法--差分方法。
先来讨论最简单的Euler 法。
为此，首先将求解区域[0,]T 离散化为若干个离散点：0110N N t t t t T -=<<<<= （i.3） 其中n t hn =，0h >称为步长。
在微积分课程中我们熟知，微商（即导数）是差商的极限。
反过来，差商就是微商的近似。
在0t t =处，在（i.1a ）中用向前差商10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到1000(,)u u hf t u -=一般地，我们有1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
下面我们用数值积分法重新导出 Euler 法以及其它几种方法。
为此，在区间1[,]n n t t +上积分常微分方程（i.1a ），得11()()(,())n n t n n t u t u t f t u t dt ++=+⎰ （i.5）用各种数值积分公式计算（i.5）中的积分，便导致各种不同的差分法。
数值分析第八章 常微分方程的数值解
Euler方法的局部截断误差： （设yn=y(xn)） Tn1 y(xn1) y(xn ) hf (xn , y(xn ))
y(xn h) y(xn ) hy(xn )
y(xn ) hy(xn )
h2 2
y( n )
y(xn ) hy(xn )
h2 2
y( n )
h2 2
y(xn ) O(h3 )
n (xn , xn1 )
其中 (h2 / 2) y(xn ) 称局部截断误差主项。
故Tn+1= O(h2)，p=1， 即Euler方法具1阶精度。
梯形方法的局部截断误差： （设yn=y(xn)）
Tn1
y(xn1 )
y(xn )
h 2
[
y
(
xn
)
y(xn1 )]
hy(xn )
h2 2
y(xn )
由Taylor展式
y(xn1) y(xn h)
yn1
y(xn ) hy(xn )
h2 2
y(xn )
hp p!
y( p) (xn )
Tn+1= O(hp+1)，若提高p，可提高精度。 但因 y f (x, y)
y(x) fx(x, y) f y(x, y) f (x, y)
……
高阶导数计算复杂，故可从另外角度考虑。
Tn1 y(xn1 ) y(xn ) h (xn , y(xn ), h)
为显式单步法在节点xn+1处的局部截断误差 。 若存在最大整数p使局部截断误差满足
Tn1 y(xn h) y(xn ) h(xn , yn , h) O(h p1)
则称显式单步法具有p阶精度或称p阶方法。 注：将Tn+1表达式各项在xn处作Taylor展开，可 得具体表达式。
常微分方程数值解法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最基础的方法之一，其基本思想是通过离散化时间点上的函数值来 逼近微分方程的解。
详细描述
欧拉方法基于微分方程的局部线性化，通过在时间点上逐步逼近微分方程的解，得到一系列离散点上 的近似值。该方法简单易行，但精度较低，适用于求解初值问题。
龙格-库塔方法
总结词
影响
数值解法的稳定性对计算结果的精度和可靠 性有重要影响。
判断方法
通过分析数值解法的迭代公式或离散化方法， 判断其是否具有稳定性和收敛性。
数值解法的收敛性
定义
数值解法的收敛性是指随着迭代次数的增加， 数值解逐渐接近于真实解的性质。
影响
数值解法的收敛性决定了计算结果的精度和 计算效率。
分类
根据收敛速度的快慢，可以分为线性收敛和 超线性收敛等。
判断方法
通过分析数值解法的迭代公式或离散化方法， 判断其是否具有收敛性。
误差分析
定义
误差分析是指对数值解法计算过程中 产生的误差进行定量分析和估计的过 程。
分类
误差可以分为舍入误差、截断误差和 初始误差等。
影响
误差分析对于提高计算精度和改进数 值解法具有重要意义。
分析方法
通过建立误差传递公式或误差估计公 式，对误差进行定量分析和估计。
生物学
生态学、生物种群动态和流行病传播 等问题可以通过常微分方程进行建模
和求解。
化学工程
化学反应动力学、化学工程流程模拟 等领域的问题可以通过常微分方程进 行描述和求解。
经济学
经济系统动态、金融市场模拟和预测 等问题可以通过常微分方程进行建模 和求解。
02 常微分方程的基本概念
常微分方程的定义
常微分方程的数值解法
常微分方程的数值解法常微分方程是研究变量的变化率与其当前状态之间的关系的数学分支。
它在物理、工程、经济等领域有着广泛的应用。
解常微分方程的精确解往往十分困难甚至不可得，因此数值解法在实际问题中起到了重要的作用。
本文将介绍常见的常微分方程的数值解法，并比较其优缺点。
1. 欧拉方法欧拉方法是最简单的数值解法之一。
它基于近似替代的思想，将微分方程中的导数用差商近似表示。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）根据微分方程的定义使用近似来计算下一个点的值。
欧拉方法的计算简单，但是由于误差累积，精度较低。
2. 改进欧拉方法为了提高欧拉方法的精度，改进欧拉方法应运而生。
改进欧拉方法通过使用两个点的斜率的平均值来计算下一个点的值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）根据微分方程的定义使用近似来计算下一个点的值。
改进欧拉方法相较于欧拉方法而言，精度更高。
3. 龙格-库塔法龙格-库塔法（Runge-Kutta）是常微分方程数值解法中最常用的方法之一。
它通过迭代逼近精确解，并在每一步中计算出多个斜率的加权平均值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）计算各阶导数的导数值。
（4）根据权重系数计算下一个点的值。
与欧拉方法和改进欧拉方法相比，龙格-库塔法的精度更高，但计算量也更大。
4. 亚当斯法亚当斯法（Adams）是一种多步法，它利用之前的解来近似下一个点的值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）通过隐式或显式的方式计算下一个点的值。
亚当斯法可以提高精度，并且比龙格-库塔法更加高效。
5. 多步法和多级法除了亚当斯法，还有其他的多步法和多级法可以用于解常微分方程。
多步法通过利用多个点的值来逼近解，从而提高精度。
而多级法则将步长进行分割，分别计算每个子问题的解，再进行组合得到整体解。
常微分方程数值解法ppt课件
若存在正的常数 L 使：
(Lipschitz)条件
|f( x ,y 1 ) f( x ,y 2 ) | L |y 1 y 2 | ( 1 .3 )
使 得 对 任 意 的 x [ a , b ] 及 y 1 ,y 2 都 成 立
则称 f (x,y) 对y 满足李普希兹条件，L 称为 Lipschitz常数.
节点 x i a i h i , 一 般 取 h i h ( ( b a ) / n ) 即 等 距
要计算出解函数 y(x) 在一系列节点
ax 0x 1x nb
处的近似值 y y(x ) i 完整版PPT课件i
16
yf(x,y) axb (1 .1 )
y(x 0) y0
(1 .2 )
对微分方程(1.1)两端从 xn到 xn1 进行积分
在大量的实际方程中出现的函数起码的连续性都 无法保证，更何况要求阶的导数
求解数值解
很多微分方程 根本求不到 问题的解析解！
重要手段。
完整版PPT课件
7
5.常微分方程数值解法的特点 常微分方程的数值解法常用来求近似解
根据提供的算法 通过计算机
数值解法得到的近似 解(含误差)是一个 离散的函数表.
便捷地实现
欧拉方法的导出把区间ab分为n个小区间步长为要计算出解函数yx在一系列节点iiyyx?iiixaihhhban?????一般取即等距节点处的近似值01naxxxb?????1iiihxx??nn等分001112yfxyaxbyxy????????对微分方程11两端从1nnxx?到进行积分11nnnnxxxxydxfxyxdx??????11nnxnnxyxyxfxyxdx?????右端积分用左矩形数值求积公式22baggxdxbagaba???????gxfxyx?令11nnnnxxnnfxyxnnyyfxyxh??????得x0x11nnnnnnyxyxhyxyhfxy??????1
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
郑州大学研究生课程（2012-2013学年第一学期）数值分析Numerical Analysis习题课第八章常微分方程数值解法待求解的问题：一阶常微分方程的初值问题/* Initial-Value Problem */:⎪⎩⎪⎨⎧=∈=0)(],[),(y a y b a x y x f dx dy 解的存在唯一性（“常微分方程”理论）：只要f (x , y ) 在[a , b ] ×R 1 上连续，且关于y 满足Lipschitz 条件，即存在与x , y 无关的常数L 使对任意定义在[a , b ] 上的y 1(x ) 和y 2(x ) 都成立，则上述IVP 存在唯一解。
1212|(,)(,)|||f x y f x y L y y −≤−一、要点回顾§8.2 欧拉(Euler)法通常取(常数),则Euler 法的计算格式h h x x i i i ==−+1⎩⎨⎧=+=+)(),(001x y y y x hf y y i i i i i =0,1,…,n ( 8.2 )§8.2 欧拉(Euler)法(1) 用差商近似导数))(,()()()()(1n n n n n n x y x hf x y x y h x y x y +=′+≈+⎩⎨⎧=+=+)(),(01a y y y x hf y y n n n n 差分方程初值问题向前Euler 方法h x y x y x y n n n )()()(1−≈′+))(,()()(1n n n n x y x f hx y x y ≈−+))(,()(n n n x y x f x y =′§8.2 欧拉(Euler)法若用向后差商近似导数，即))(,()()(111++++≈n n n n x y x hf x y x y ⎩⎨⎧=+=+++)(),(0111a y y y x hf y y n n n n hx y x y x y n n n )()()(11−≈′++向后Euler 方法))(,()()(111+++≈−n n n n x y x f h x y x y ))(,()(111+++=′n n n x y x f x y§8.2 欧拉(Euler)法（2）用数值积分方法∫+=−+1))(,()()(1n n x x n n dxx y x f x y x y ∫∫++=′11))(,()(n n n n x x x x dx x y x f dx x y 1111(,)(,()), (,)(,()),()()n n n n n n n n f x y f x y x f x y f x y x y x y y x y ++++≈≈≈≈分别用左矩形和右矩形公式，即代替上式右端的积分,并注意 ,分别得到1111(,)(,)n n n n n n n n y y h f x y y y h f x y ++++=+=+，。
向前欧拉公式和向后欧拉公式：§8.2 欧拉(Euler)法若对积分用梯形公式，则得[]))(,())(,(2)()(111++++≈−n n n n n n x y x f x y x f h x y x y []⎪⎩⎪⎨⎧=++≈+++)(),(),(20111a y y y x f y x f h y y n n n n n n 梯形欧拉公式§8.2 欧拉(Euler)法欧拉方法的收敛性11212 ()Taylor , (,),()()()()()2()(,())() (1)2n n n n n n n n n n n n n n y x x x x h y x y x h y x hy x y h y x hf x y x y ξξξ+++∃∈′′′=+=++′′=++将在点展开%1(),(8.2) () (,()) (2)n n n nn y x y y x h f x y x ++ 假定已知准确值 利用欧拉公式，定义§8.2 欧拉(Euler)法局部截断误差%()11112() ()()(,()) (). 2n n n n n n n n T y x y y x y x hf x y x h y ξ++++=−=−+′′=2221 max |()|,|||()|(). 22a x bn n M y x h h T y M O h ξ≤≤+′′=′′=≤=令则称为局部截断误差§8.2 欧拉(Euler)法欧拉方法的收敛性定义若给定方法的局部截断误差满足则称该方法是P 阶的，或称为具有P 阶精度。
11||(),p n T O h ++=§8.2 欧拉(Euler)法整体截断误差111().n n n e y x y +++=−记112121112 ,,(), (),, () ,(),,,.n n n n n n y y y y y x y x y x y x e y y y +++L L L 因为计算 时用到的,,是 的近似值每步产生的误差会累积到计算的误差中因此 与 ,, 都有关称整体截断误差为%%%1111111111|||()||()||| ||||. (3)n n n n n n n n n n e y x y y x y y y T y y ++++++++++=−≤−+−≤+−§8.2 欧拉(Euler)法欧拉方法的收敛性%%()1111,(8.2)|||()(,())(,)| |() ||(,())(,)| |() ||() | () (1)n n n n n n n n n n n n n n n n n n n n y y y y y x h f x y x y h f x y y x y h f x y x f x y y x y hL y x y hL ++++−=+−+≤−+−≤−+−≤+对应用欧拉公式得李普希兹条件||. (4)n e由此知，当max ||,(4)(3)k kT T =记 将代入得[]1121210112||(1)||(1)(1)|| (1)(1)|| (1)(1)(1) (1)||(1)1(1)1()(1)1 (n n n n nn n n e T hL e T hL T hL e T hL T hL e T hL T hL T hL T hL e hL hL T O h hL hLO +−−+++≤++≤++++≤++++≤≤++++++++++−+−≤=+−≤L L ).h 10,||0, n h e +→→欧拉公式是一阶有收敛的。
§8.2 欧拉(Euler)法向后欧拉公式#隐式欧拉法或向后欧拉法#/* implicit Euler method or backward Euler method*/11()()()n n n y x y x y x h++−′≈x n +1点向后差商近似导数111111()()()()()(,)n n n n n n n n n n y x y x hy x y x y y x y y h f x y ++++++′≈+↑≈≈=+代入隐式或后退欧拉公式§8.2 欧拉(Euler)法向后欧拉公式由于未知数y n +1同时出现在等式的两边，故称为隐式/* implicit */欧拉公式，而前者称为显式/* explicit */ 欧拉公式。
隐式公式不能直接求解，一般需要用Euler 显式公式得到初值，然后用Euler 隐式公式迭代求解。
因此隐式公式较显式公式计算复杂，但稳定性好（后面分析）。
隐式欧拉公式中的未知数y n +1可通过以下迭代法求解：0)1(1)()111(,)(,)n n n n k k n n n n y y h f x y yy h f x y +++++⎧=+⎨=+⎩（§8.2 欧拉(Euler)法向后欧拉公式1(,) 0,1,...n n n n y y h f x y n +=+=比较欧拉显式公式和隐式公式及其局部截断误差231112()()()h n n n n T y x y y x O h +++′′=−=+%显式公式111(,)n n n n y y h f x y +++=+隐式公式231112()()()h n n n n T y x y y x O h +++′′=−=−+%若将这两种方法进行算术平均，即可消除误差的主要部分/*leading term*/而获得更高的精度,称为梯形法#梯形公式/* trapezoid formula */—显、隐式两种算法的平均111[(,)(,)]2n n n n n n hy y f x y f x y +++=++R y xyO h+++=−=例8.2.3对初值问题⎩⎨⎧==+′1)0(0y y y 证明用梯形公式求得的近似解为nn h h y ⎟⎠⎞⎜⎝⎛+2−2=并证明当步长h →0时,y n 收敛于精确解证明: 解初值问题的梯形公式为xe−)],(),([1+1+1++2+=n n n n n n y x f y x f hy y y y x f −=),(∵][211++−−+=n n n n y y hy y ∴整理成显式n n y h h y ⎟⎠⎞⎜⎝⎛+2−2=1+反复迭代,得到1+2−31−21+⎟⎠⎞⎜⎝⎛+2−2==⎟⎠⎞⎜⎝⎛+2−2=⎟⎠⎞⎜⎝⎛+2−2=⎟⎠⎞⎜⎝⎛+2−2=y h h y h h y h h y h h y n n n n n ...10=y 22nxnh y e h −−⎛⎞=→⎜⎟+⎝⎠∵∴梯形公式单步好隐1阶欧拉隐式公式单步差显1阶欧拉显式公式步数稳定性显隐精度局部截断误差公式()3(3)3nhy x()2(2)2nhy x()2(2)2nhy x−欧拉法小结单步好隐2阶§8.3 改进欧拉(Euler)方法先用欧拉公式(8.2)求出一个初步的近似值,称为预测值, 它的精度不高, 再用梯形公式对它校正一次,即迭代一次,求得y n+1,称为校正值, 这种预测-校正方法称为改进的欧拉公式：1n y +[]⎪⎩⎪⎨⎧++=+=++++校正预测 ),(),(2),(1111n n n n n n n n n n y x f y x f hy y y x hf y y 称为Euler 公式与梯形公式的预测—校正系统。
predictor corrector),(1iiiiyxhfyy+=+[]),(),(2111+++iiiiyxfyxf[]),(),(211 1++ +++=iiiiiiyxfyxfhyy§8.3 改进欧拉(Euler)方法可以证明,改进的欧拉公式的精度为二阶。