常微分方程数值解法共7页文档
常微分方程数值解法
用分段的折线逼近函数,此为 “折线法”而非“切线法”， 除第一个点是曲线上的切线，
其它都不是。
2、Euler方法的误差估计
1)局部截断误差。 在一步中产生的误差而非累积误差：
~
T x y y
n1
n1
n1
其中
~
y
是当
y
n
y(
x
)
n
（精确解！）时
n1
由Euler法求出的值，即y 无误差！ n
y y x , y x y 则得： h f
f ,
n1
n2
n
n
n1
n1
同样与Euler法结合，形成迭代算法，对n 0，1，2，
y y x , y 0 hf
n1
n
n
n
y y x , y x y
k 1
推出总体误差与步长的关系。
由微分方程解的存在唯一性,自然假定 （ f x，y）
充分光滑，或满足 Lipschitz条件：
f
x,ny源自xn
f
x
,
n
y
n
L
yxn
y n
第 n 步 的 总 体 截 断 误 差 记 为
en y
xn
y n
则 对 n 1 步：
e x y x y y y T y y ~ ~
用yn1, yn代替y(xn1), y(xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f
(x,
y)dx
hf
(xn ,
yn )
数值分析第九章常微分方程数值解法
松弛法
通过迭代更新函数值并逐步放松约束 条件来逼近解，适用于刚性和非刚性 问题。
利用线性组合迭代函数值来逼近解， 具有更高的收敛速度和稳定性。
03
数值解法的稳定性分析
数值解法的稳定性定义
数值解法的稳定性是指当微分方程的初值有微小的扰动时， 其数值解的近似值的变化情况。如果数值解在微小扰动下变 化较小，则称该数值方法是稳定的。
更高的精度和稳定性。
数值逼近法
泰勒级数法
将微分方程的解展开为泰勒级数，通过截断级数来逼 近解。
多项式逼近法
利用多项式来逼近微分方程的解，通过选取合适的基 函数和系数来提高逼近精度。
样条插值法
利用样条函数来逼近微分方程的解，具有更好的光滑 性和连续性。
迭代法
雅可比迭代法
通过迭代更新函数值来逼近微分方程 的解，具有简单易行的优点。
初值和边界条件的处理
根据实际问题，合理设定初值和边界 条件，以获得更准确的数值解。
收敛性和误差分析
对数值解进行收敛性和误差分析，评 估解的精度和稳定性。
数值解法的应用案例分析
人口增长模型
通过数值解法求解人口增长模型，预测未来人口数量，为政策制 定提供依据。
化学反应动力学
利用数值解法研究化学反应的动力学过程，模拟反应过程和结果。
数值分析第九章常微分方 程数值解法
• 引言 • 常微分方程数值解法的基本思想 • 数值解法的稳定性分析 • 数值解法的收敛性和误差分析 • 数值解法的实现和应用案例
01
引言
常微分方程的应用背景
自然科学
描述物理、化学、生物等自然 现象的变化规律。
工程领域
控制系统设计、航天器轨道计 算等。
第8章常微分方程数值解法
的解为
y ( x) e
x2
x 0
e dt
t2
但要计算它的值，还需要用数值积分的方法。如果要 对许多个 x 值计算解 y(x) 的近似值，那么工作量非常大。况 且实际计算不一定要求解析表达式，而是只需求在某些点 上满足精度的解的近似值或解的近似表达式就可以了。
由于高阶常微分方程可以转化为一阶常微分方程组，因 此，为了不失一般性，本章主要介绍一类一阶常微分方程初 值问题
的解来近似微分方程初值问题(8.2)的解，其 中 h (b- a) / 2 ，式(8.3)也称为欧拉公式。
欧拉法的几何意义是用一条自点 ( x0 , y0 ) 出发的 折线去逼近积分曲线 y f (x) ，如图8.1所示。 因此，这种方法又称为折线法。显然，欧拉法 简单地取折线的端点作为数值解，精度非常差。
float euler(float x0,float xn,float y0,int N) { float x,y,h; int i; x=x0; y=y0; h=(xn-x0)/(float)N; /* 计算步长 */ for(i=1;i<=N;i++) /* 欧拉公式 */ { y=y+h*func(x,y); x=x0+i*h; } return(y); }
8.4 龙格—库塔（Runge-Kutta）法 8.4.1 龙格—库塔法的基本思想
在欧拉法 yi 1 yi h f ( xi , yi ) (i 0,1,) 中，用解函数 y f (x) 在 点 x i 处的斜率 f ( xi , y i ) 计算从 yi 到 y i 1 的增量，y i 1 的表达式 与 y( xi 1 ) 的Taylor展开式的前二项相等，使方法只有一阶精度。 改进的欧拉法用两个点 x i ，x i 1 处的斜率 f ( xi , y i )、f ( xi 1 , yi 1 ) 的平均值计算增量，使方法具有二阶精度，即 y i 1 的表达式 与 y( xi 1 ) 的Taylor展开式的前三项相等。 由此龙格和库塔提出了一种间接地运用Taylor公式的方法， y (x) 即利用 在若干个待定点上的函数值和导数值做出线性组 合式，选取适当系数使这个组合式进行Taylor展开后与 y( xi 1 ) 的Taylor展开式有较多的项达到一致，从而得出较高阶的数 值公式，这就是龙格—库塔法的基本思想。
数值分析常微分方程数值解法
第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
第十章.常微分方程数值解法
n+1 n k +1 n+1 n
(
n
n
k
n+1
n+1
)
k = 0, ,2,L 1
（二）向后Euler 法的稳定性 对y = λy 用 后 向 Euler法： = y + hλ y y 误 公 ： 差 式 ρ = ρ + hλ ρ ρ = 1 = 1 则 1− hλ ρ [(1−λh)(1−λh)]
( 0)
xn
(
(
))
(
)
(
)
梯形公式的稳定性
梯形公式比Euler 方法的局部与总体误差均高一阶， = 0(h )， e
2 N
但每次迭代均多算一次函数值—提高精度的计算代价
ρ
n+ 1
= ρn +
λh
2
(ρ + ρ )
n n+ 1 1 2
常微分方程的数值解法
第六章 常微分方程的数值解法 §6.0 引言§6.1 算法构造的主要途径§6.2 Runge-Kutta Method 算法§6.3 线性多步法§6.4 线性多步法的一般形式§6.5 算法的稳定性、收敛性§6.0 引 言1． 主要考虑如下的一阶常微分方程初值问题的求解:()()00,dy f x y dx y x y ⎧=⎪⎪⎨=⎪⎪⎩ 微分方程的解就是求一个函数y=y(x)，使得该函数满足微分方程并且符合初值条件。
2. 例如微分方程:xy '-2y=4x ；初始条件: y(1)=-3。
于是可得一阶常微分方程的初始问题24(1)3y y x y ⎧'=+⎪⎨⎪=-⎩。
显然函数y(x)=x 2-4x 满足以上条件,因而是该初始问题的微分方程的解。
3. 但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示。
因此，只能依赖于数值方法去获得微分方程的数值解。
4． 微分方程的数值解：设微分方程问题的解y(x)的存在区间是[a,b ]，初始点x 0=a ，将[a,b ]进行划分得一系列节点x 0 , x 1 ,...,x n ，其中a= x 0< x 1<…< x n =b 。
y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点x k 的近似值y(x k )，即y≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。
如图所示:§6.1 算法构造的主要途径x 0 x 1 x 2 ...1 欧拉公式1.1 构造的思想：利用差商代替一阶导数，即010()()x x y x y x dy dx h =-≈，则 1000()()(,)y x y x f x y h -≈。
第7章 常微分方程数值解法
代入(6―3)式得
h yi 1 yi [ f ( xi , yi ) f ( xi 1 , yi 1 )] 2 i 0,1, 2, , n 1
(6―5)
这样得到的点列仍为一折线，只是用平均斜率 来代替原来一点处的斜率。式(6―5)称为改进的欧拉 公式。
不难发现，欧拉公式(6―3)是关于yi+1 的显式，只
h y xi 1 yi 1 f xi 1 , y xi 1 f xi 1 , yi 1 2 (6―15) h 3 '' f 12
因此
hL h3 y ( xi 1 ) yi 1 y ( xi 1 ) yi 1 f ( ) 2 12 h3 (1 q) y ( xi 1 ) yi 1 f ( ) 12 y ( xi 1 ) yi 1 O ( h 3 )
c 并取 yi 1 yi(1)
(6―7)
虽然式(6―7)仅迭代一次，但因进行了预先估计，
故精度却有较大的提高。 在实际计算时，还常常将式(6―7)写成下列形式：
k1 f ( xi , yi ) k f ( x h, y hk ) i i 1 2 h yi 1 yi 2 (k1 k2 ) i 0,1, 2,
在进行误差分析时，我们假设yi=y(xi)，考虑用
yi+1 代替y(x
i+1)而产生截断误差，确定欧拉公式和改
进的欧拉公式的精确度。 设初值问题(6―1)的准确解为y=y(x)，则利用泰 勒公式
y ( xi 1 ) y ( xi h ) h2 y ( xi ) hy ( xi ) y ( xi ) 2! h3 y ( xi ) 3!
常微分方程初值问题数值解法
0.4 1.3582 1.3416 0.9 1.7178 1.6733
0.5 1.4351 1.4142 1.0 1.7848 1.7321
7
初值问题（2.2）有解 y ，1按2这x 个解析式子
算出的准确值 y(x同n )近似值 一y起n 列在表9-1中，两者 相比较可以看出欧拉方法的精度很差.
17
所以，局部截断误差可理解为用方法（2.10）计算一步的 误差，也即公式（2.10）中用准确解y(x代) 替数值解产生
的公式误差.
根据定义，显然欧拉法的局部截断误差
Tn1 y( xn1) y( xn ) hf ( xn , y( xn ))
y(xn h) y(xn ) hy(xn )
y(2) n1
yn
hf
( xn1,
y (1) n1
).
11
如此反复进行，得
y (k 1) n1
yn
hf
( xn1,
y(k) n1
),
(k 0,1, ).
（2.6）
由于 f (x,对y) 满足y 利普希茨条件（1.3）. 由（2.6）减 （2.5）得
y (k 1) n 1
yn1
h
f
( xn1,
y(k) n 1
积分曲线上一点 (x的, y切)线斜率等于函数 值.
的f (x, y)
如果按函数 f (在x, y) 平x面y上建立一个方向场，那 么，积分曲线上每一点的切线方向均与方向场在该点的方 向相一致.
基于上述几何解释，从初始点 P0 (x出0 ,发y0，) 先依 方向场在该点的方向推进到 x 上x1一点 ，P然1 后再从 P1 依方向场的方向推进到 x 上x2一点 ，循P2此前进做出
第5章常微分方程数值解法
2hyxn
2h3 3!
y
yxn1
yxn1 2hyxn
h3 3
y ②
将①、②两式相减：
y
xn1
h3 yn1 3
y
——两步法局部截断误差
18
第19页/共24页
2024年8月7日
（2）梯形公式
yn1
yn
h 2
f
xn , yn
f
xn1 , yn1
yxn
hf
2
xn ,
yn1 yxn1 2hyxn ①
2024年8月7日
17
第18页/共24页
2024年8月7日
将函数用泰勒级数展开：（ h 较小， 相差不大）
yxn1
yxn
h 1!
yxn
h2 2!
yxn
h3 3!
y
yxn1
yxn
h
1!
yxn
h2
2!
yxn
h3
3!
y
yxn1
yxn1
yxn
f
xn1 , yxn1
yn1
yxn
h 2
y
xn
yxn1 ①
将函数用泰勒级数展开：
yxn1
yxn
h 1!
yxn
h2 2!
yxn
h3 3!
y ②
yxn1
yxn
h 1!
yxn
h2 2!
y ③
（ h 较小， 相差不大）
19
第20页/共24页
2024年8月7日
①、②两式相减，并代入③式：
（图示表示梯形法计算结果）
常微分方程的数值解法全文
第8章常微分方程的数值解法8.4单步法的收敛性与稳定性8.4.1相容性与收敛性上面所介绍的方法都是用离散化的方法，将微分方程初值问题化为差分方程初值问题求解的.这些转化是否合理？即当h →∞时，差分方程是否能无限逼近微分方程，差分方程的解n y 是否能无限逼近微分方程初值问题的准确解()n y x ，这就是相容性与收敛性问题.用单步法(8.3.14)求解初值问题(8.1.1)，即用差分方程初值问题100(,,)()n n n n y y h x y h y x y ϕ+=+⎧⎨=⎩(8.4.1)的解作为问题(8.1.1)的近似解，如果近似是合理的，则应有()()(,(),)0 (0)y x h y x x y x h h hϕ+--→→(8.4.2)其中()y x 为问题(8.1.1)的精确解.因为0()()lim ()(,)h y x h y x y x f x y h→+-'==故由(8.4.2)得lim (,,)(,)h x y h f x y ϕ→=如果增量函数(,(),)x y x h ϕ关于h 连续，则有(,,0)(,)x y f x y ϕ=(8.4.3)定义8.3如果单步法的增量函数(,,)x y h ϕ满足条件(8.4.3)，则称单步法(8.3.14)与初值问题(8.1.1)相容.通常称(8.4.3)为单步法的相容条件.满足相容条件(8.4.3)是可以用单步法求解初值问题(8.1.1)的必要条件.容易验证欧拉法和改进欧拉法均满足相容性条件.一般地，如果单步法有p 阶精度（1p ≥），则其局部截断误差为[]1()()(,(),)()p y x h y x h x y x h O h ϕ++-+=上式两端同除以h ，得()()(,,)()p y x h y x x y h O h hϕ+--=令0h →，如果(,(),)x y x h ϕ连续，则有()(,,0)0y x x y ϕ'-=所以1p ≥的单步法均与问题(8.1.1)相容.由此即得各阶龙格－库塔法与初值问题(8.1.1)相容.定义8.4一种数值方法称为是收敛的，如果对于任意初值0y 及任意固定的(,]x a b ∈，都有lim () ()n h y y x x a nh →==+其中()y x 为初值问题(8.1.1)的精确解.如果我们取消局部化假定，使用某单步法公式，从0x 出发，一步一步地推算到1n x +处的近似值1n y +.若不计各步的舍入误差，而每一步都有局部截断误差，这些局部截断误差的积累就是整体截断误差.定义8.5称111()n n n e y x y +++=-为某数值方法的整体截断误差.其中()y x 为初值问题(8.1.1)的精确解，1n y +为不计舍入误差时用某数值方法从0x 开始，逐步得到的在1n x +处的近似值（不考虑舍入误差的情况下，局部截断误差的积累）.定理8.1设单步法(8.3.14)具有p 阶精度，其增量函数(,,)x y h ϕ关于y 满足利普希茨条件，问题(8.1.1)的初值是精确的，即00()y x y =，则单步法的整体截断误差为111()()p n n n e y x y O h +++=-=证明由已知，(,,)x y h ϕ关于y 满足利普希茨条件，故存在0L >，使得对任意的12,y y 及[,]x a b ∈，00h h <≤，都有1212(,,)(,,)x y h x y h L y y ϕϕ-≤-记1()(,(),)n n n n y y x h x y x h ϕ+=+，因为单步法具有p 阶精度，故存在0M >，使得1111()p n n n R y x y Mh ++++=-≤从而有111111111()()()(,(),)(,,)()(,(),)(,,)n n n n n n n p n n n n n n p n n n n n n e y x y y x y y y Mh y x h x y x h y h x y h Mh y x y h x y x h x y h ϕϕϕϕ+++++++++=-≤-+-≤++--≤+-+-1(1)p nMh hL e +≤++反复递推得11111101110(1)(1)1(1)(1)(1)(1)1(1)p p n n n p n n p n e Mh hL Mh hL e hL hL Mh hL e hL Mh hL e hL+++-+++++⎡⎤≤++++⎣⎦⎡⎤≤+++++++⎣⎦+-≤++因为00()y x y =，即00e =，又(1)n h b a +≤-，于是ln(1)1()(1)(1)b a b a hL n L b a h h hL hL e e --++-+≤+=≤所以()11()p L b a p n M e h e O h L -+⎡⎤≤-=⎣⎦推论设单步法具有p (1p ≥)阶精度，增量函数(,,)x y h ϕ在区域G ：, , 0a x b y h h ≤≤-∞<<+∞≤≤上连续，且关于y 满足利普希茨条件，则单步法是收敛的.当(,)f x y 在区域:,D a x b y ≤≤-∞<<+∞上连续，且关于y 满足利普希茨条件时，改进欧拉法，各阶龙格－库塔法的增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，因而它们都是收敛的.关于单步法收敛的一般结果是：定理8.2设增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，则单步法收敛的充分必要条件是相容性条件(8.4.3).8.4.2稳定性稳定性与收敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长0h →时，方法的整体截断误差是否趋于零的问题.而稳定性则是讨论舍入误差的积累能否对计算结果有严重影响的问题.定义8.6若一种数值方法在节点值n y 上有一个大小为δ的扰动，于以后各节点()m y m n >上产生的偏差均不超过δ，则称该方法是稳定的.我们以欧拉法为例进行讨论.假设由于舍入误差，实际得到的不是n y 而是n n n y y δ=+，其中n δ是误差.由此再计算一步，得到1(,)n n n n y y hf x y +=+把它与不考虑舍入误差的欧拉公式相减，并记111n n n y y δ+++=-，就有[]1(,)(,)1(,)n n n n n n y n nh f x y f x y hf x δδηδ+⎡⎤=+-=+⎣⎦其中y f f y∂=∂.如果满足条件1(,)1y n hf x η+≤，(8.4.4)则从n y 到1n y +的计算，误差是不增的，可以认为计算是稳定的.如果条件(8.4.4)不满足，则每步误差将增大.当0y f >时，显然条件(8.4.4)不可能满足，我们认为问题本身具有先天的不稳定性.当0y f <时，为了满足稳定性要求(8.4.4)，有时h 要很小.一般的，稳定性与方法有关，也与步长h 的大小有关，当然也与方程中的(,)f x y 有关.为简单起见，通常只考虑数值方法用于求解模型方程的稳定性，模型方程为y y λ'=(8.4.5)其中λ为复数.一般的方程可以通过局部线性化转化为模型方程，例如在(,)x y 的邻域(,)(,)(,)()(,)()x y y f x y f x y f x y x x f x y y y '==+-+-+略去高阶项，再作变量替换就得到u u λ'=的形式.对于模型方程(8.4.5)，若Re 0λ>，类似以上分析，可以认为方程是不稳定的.所以我们只考虑Re 0λ<的情形，这时不同的数值方法可能是数值稳定的或者是数值不稳定的.当一个单步法用于试验方程y y λ'=，从n y 计算一步得到1()n n y E h y λ+=(8.4.6)其中()E h λ依赖于所选的方法.因为通过点(,)n n x y 试验方程的解曲线（它满足,()n n y y y x y λ'==）为[]exp ()n n y y x x λ=-，而一个p 阶单步法的局部截断误差在()n n y x y =时有1111()()p n n n T y x y O h ++++=-=，所以有1exp()()()p n n y h E h y O h λλ+-=(8.4.7)这样可以看出()E h λ是h e λ的一个近似值.由(8.4.6)可以看到，若n y 计算中有误差ε，则计算1n y +时将产生误差()E h λε，所以有下面定义.定义8.7如果(8.4.6)式中，()1E h λ<，则称单步法(8.3.14)是绝对稳定的.在复平面上复变量h λ满足()1E h λ<的区域，称为方法(8.3.14)的绝对稳定区域，它与实轴的交称为绝对稳定区间.在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致.事实上，()1E h λ=时也可以认为误差不增长.(1)欧拉法的稳定性欧拉法用于模型方程(8.4.5)，得1(1)n n y h y λ+=+，所以有()1E h h λλ=+.所以绝对稳定条件是11h λ+<，它的绝对稳定区域是h λ复平面上以(1,0)-为中心的单位圆，见图8.3.而λ为实数时，绝对稳定区间是(2,0)-.Im()h λRe()h λ2-1-O 图8.3欧拉法的绝对稳定区域(2)梯形公式的稳定性对模型方程，梯形公式的具体表达式为11()2n n n n h y y y y λλ++=++，即11212n nh y y h λλ++=-，所以梯形公式的绝对稳定区域为12112h h λλ+<-.化简得Re()0h λ<，因此梯形公式的绝对稳定区域为h λ平面的左半平面，见图8.4.特别地，当λ为负实数时，对任意的0h >，梯形公式都是稳定的.Im()h λRe()h λO 图8.4梯形公式的绝对稳定区域(3)龙格－库塔法的稳定性与前面的讨论相仿，将龙格－库塔法用于模型方程(8.4.5)，可得二、三、四阶龙格－库塔法的绝对稳定区域分别为211()12h h λλ++<23111()()126h h h λλλ+++<2341111()()()12624h h h h λλλλ++++<当λ为实数时，二、三、四阶显式龙格－库塔法的绝对稳定区域分别为20h λ-<<、2.510h λ-<<、 2.780h λ-<<.例8.5设有初值问题21010101(0)0xy y x x y ⎧'=-≤≤⎪+⎨⎪=⎩用四阶经典龙格－库塔公式求解时，从绝对稳定性考虑，对步长h 有何限制？解对于所给的微分方程有2100,(010)1f x x y xλ∂==-<≤≤∂+在区间[0,10]上，有201010max ||max51t x x λ<<==+由于四阶经典龙格－库塔公式的绝对稳定区间为 2.7850h λ-<<，则步长h 应满足00.557h <<.。
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解y′=f(x,y),x∈[x0,b]y(x0)=y0.定理1（利普希茨条件）若存在正数L，使得对任意，y1，y2，有|f(x,y1)−f(x,y2)|≤L|(y1−y2)|定理2（解存在性）①若函数f在方区域x∈[a,b],y∈R连续，②函数f关于y 满足利普希茨条件，则对任意x∈[a,b]，常微分方程存在唯一的连续可微数值解.两类问题：①单步法---计算下一个点的值yn+1只需要用到前面一个点的值yn②多步法---计算下一个点的值yn+1需要用到前面l个点的值yl1、欧拉法---下一个点的计算值等于前一个点的计算值加上步长乘以前一个点的函数值•具体过程一些批注：显式欧拉方程指下一步要计算的值，不在迭代方程中；隐式欧拉方程指下一步要计算的值，在迭代方程中。
怎么计算隐式欧拉方程----要借助显示欧拉迭代计算---一般用迭代法-----迭代---将微分方程在区间[xn,xn+1]进行积分，然后函数f进行近似，即可得到迭代方程-----迭代方程收敛性？由函数关于y满足利普希茨条件，可以推出迭代公式收敛。
•局部截断误差：假设前n步误差为0，我们计算第n+1步的误差，将次误差称为局部截断误差，且局部误差为O(hp+1)•p阶精度：由理论证明：若局部误差阶的时间复杂度为O(hp+1)，则整体误差阶为O(hp)我们称公式精度为p。
•显示欧拉法与隐式欧拉法•梯形方法----将显式欧拉迭代方程与隐式欧拉迭代方程做一下加权平均，构造的计算公式.•改进的欧拉方法---思想：因为梯形公式是隐式公式，将显式欧拉公式对下一步的计算值进行预估，用梯形公式对下一步的计算值进行校正.2、龙格-库塔方法思想：根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以前一个点的斜率；而这个斜率用该区间上的多个点的斜率的算数平均来逼近。
注意：怎么计算任意斜率Ki？第i个点的斜率Ki有微分方程可以算出f′=f(xn,yn)所以要算的f(xn,yn)值，由欧拉法即可算出, yn+1=yn+hf′•2阶-龙格-库塔方法----类似改进的欧拉法根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
第6章 常微分方程数值解法-文档资料
for i=1:5 y=(1-h)*y+h*x x=x+h; end
科大研究生学位课程
数值分析
6.1.2 欧拉格式的改进
对微分方程 y'd yf(x,y) , x [a,b] dx
所以，有格式为：
yi 1yih 2[f(xi,yi)f(xi 1,yi 1)]
上式称为梯形格式。
类似，可以算出梯形格式的误差估计式：
i1 O(h3)
2阶的方法
梯形法是二阶、隐式单步的方法，要用迭代法求解。怎么求？
科大研究生学位课程
数值分析
改进欧拉格式 /* modified Euler’s Formula */
科大研究生学位课程
数值分析
微分方程的数值方法。求的是在一系列离散点列上，未知 函数y(x)在这些点上的值的近似。 为了考察数值方法提供的数值解，是否有实用价值，需要 知道如下问题：
① 步长充分小时，所得到的数值解能否逼近问题得真解； 即收敛性问题 ② 误差估计，主要考虑局部截断误差。 ③ 产生的舍入误差，在以后得各步计算中，是否会无限
计算定积分，则： d xi1 y xi1 f(x,y)dx
xi
xi
y(xi1)y(xi)xxii1 f(x,y)dx
用左矩形公式 x x ii 1f(x ,y (x )d ) x f(x i,y (x i)()x i 1 x i)
y ( x i 1 ) y i f( x i,y ( x i)x ) i 1 ( x i)
数值分析
例6.1：考虑初值问题：
计算方法 第6章 常微分方程数值解
已知Euler格式 yn1 yn hf ( xn , yn )
h2 y( xn1 ) yn1 2 y''( xn )
即Euler格式具有一阶精度
如果令
y( xn1 ) y( xn1 ) 2h
y'( xn )
f ( xn , yn )
并假定 y( xn1 ) yn1, y( xn ) yn
常微分方程数值解
常微分方程的数值解法
§1 引 言 §2 欧拉方法 §3 龙格-库塔方法
2
§1 引 言
在工程和科学技术的实际问题中，常需要解常微 分方程。但常微分方程组中往往只有少数较简单和典 型的常微分方程（例如线性常系数常微分方程等）可 求出其解析解。对于变系数常微分方程的解析求解就 比较困难，而一般的非线性常微分方程就更不用说了。 在大多数情况下，常微分方程只能用近似法求解。这 种近似解法可分为两大类：一类是近似解析法，如级 数解法、逐次逼近法等；另一类则是数值解法，它给 出方程在一些离散点上的近似解。
yn
2
xn yn
令 h 0.1 将 x0 0, y0 1 代入Euler格式
步进计算结果见P106表5.1
第五章：常微分方程数值解
Euler值
y 1 2x
第五章：常微分方程数值解
Euler格式的误差分析
pn1
事实上Euler格式的每一步都存在误差，为了方便讨论y算( x)
d2x
dt 2 x(t
0
c
m )
x x
0
0 (t
t) 0
x(t ) x
0
0
5
第七章 常微分方程数值解法.
y f ( x, y), x [a,b]
y(
x0
)
y0
23
7.1 欧拉法和改进的欧拉法
欧拉公式
yi1 yi y0 y( x0 )
h
f
(xi ,
yi )
,
i
0,1,2,
改进的欧拉公式
yi
1
26
引言：
公式构造思想：从泰勒公式出发，寻找更高阶的 数值公式。
例如，泰勒公式计算到二阶可得
y(x + h) = y(x) + yⅱ(x)h + 1 y ?(x)h2 + O(h3) 2!
因
ìïïïíïïïî
y¢(x) = yⅱ(x) =
f (x, y(x)) df (x, y(x))
dx
=
fx (x, y(x)) +
可证明预测－校正公式的截断误差也为 O(h3)。
18
7.1.2 改进的欧拉法及预测-校正公式
例 取步长h=0.2,用改进的欧拉法的预测－校正公
式求解初值问题的数值解y1 , y2 .
ìïïíïïî
y¢= x + y(0) = 1
y
解
f (x, y) = x + y, x0 = 0, y0 = 1
预测-校正公式具体是
y( p) 2
=
0.2x1 +
1.2 y1
=
0.2?
0.2
1.2? 1.24
1.528
y2 = y1 + 0.1[(x1 + y1) + (x2 + y2( p) )] = 1.24 + 0.1(0.2 + 1.24 + 0.4 + 1.528) = 1.5768
第六章常微分方程的数值解法
第六章常微分方程的数值解法第六章常微分方程的数值解法在自然科学研究和工程技术领域中，常常会遇到常微分方程的求解问题。
传统的数学分析方法仅能给出一些简单的、常系数的、经典的线性方程的解析表达式，不能处理复杂的、变系数的、非线性方程，对于这些方面的问题，只能求诸于近似解法和数值解法。
而且在许多实际问题中，确确实实并不总是需要精确的解析解，往往只需获得近似的解或者解在若干个点上的数值即可。
在高等数学课程中介绍过的级数解法和逐步逼近法，能够给出解的近似表达式，这一类方法称为近似解法。
还有一类方法是通过计算机来求解微分方程的数值解，给出解在一些离散点上的近似值，这一类方法称作为数值方法。
本章主要介绍常微分方程初值问题的数值解法，包括Euler 方法、Runge-Kutta 方法、线性多步法以及微分方程组与高阶微分方程的数值解法。
同时，对于求解常微分方程的边值问题中比较常用的打靶法与有限差分法作了一个简单的介绍。
§1 基本概念1.1 常微分方程初值问题的一般提法常微分方程初值问题的一般提法是求解满足如下条件的函数，，b x a x y ≤≤)(=<<=α)(),(a y bx a y x f dxdy， (1.1) 其中),(y x f 是已知函数，α是给定的数值。
通常假定上面所给出的函数),(y x f 在给定的区域},),{(+∞<≤≤=yb x a y x D 上面满足如下条件：(1) 函数),(y x f 在区域D 上面连续；(2) 函数),(y x f 在区域D 上关于变量y 满足Lipschitz(李普希茨)条件：212121,),(),(y y b x a y y L y x f y x f ?≤≤?≤?，， (1.2)其中常数L 称为Lipschitz(李普希茨)常数。
由常微分方程的基本理论可以知道，假如(1.1)中的),(y x f 满足上面两个条件，则常微分方程初值问题(1.1)对于任意给定的初始值α都存在着唯一的解，，b x a x y ≤≤)(并且该唯一解在区间[a,b]上是连续可微的。
常微分方程数值解法
第八章 常微分方程数值解法教学目的 1. 掌握解常微分方程的单步法：Euler 方法、Taylor 方法和Runge-Kutta 方法；2. 掌握解常微分方程的多步法：Adams 步法、Simpson 方法和Milne 方法等；3. 了解单步法的收敛性、相容性与稳定性；多步法的稳定性。
教学重点及难点 重点是解常微分方程的单步法：Euler 方法、Taylor 方法和Runge-Kutta 方法和解常微分方程的多步法：Adams 步法、Simpson 方法和Milne 方法等；难点是理解单步法的收敛性、相容性与稳定性及多步法的稳定性。
教学时数 20学时 教学过程§1基本概念1．1常微分方程初值问题的一般提法常微分方程初值问题的一般提法是求函数b x a x y ≤≤),(，满足⎪⎩⎪⎨⎧=<<=)2.1()()1.1(),,(αa yb x a y x f dx dy其中),(y x f 是已知函数，α是已知值。
假设),(y x f 在区域},),{(+∞<≤≤=y b x a y x D 上满足条件： （1）),(y x f 在D 上连续； （2）),(y x f 在D 上关于变量y 满足Lipschitz 条件：2121),(),(y y L y x f y x f -≤-，21,,y y b x a ∀≤≤ （1.3）其中常数L 称为Lipschitz 常数。
我们简称条件（1）、（2）的基本条件。
由常微分方程的基本理论，我们有：定理1 当),(y x f 在D 上满足基本条件时，一阶常微分方程初值问题（1.1）、（1.2）对任意给定α存在唯一解)(x y 在],[b a 上连续可微。
定义1 方程（1.1）、（1.2）的解)(x y 称为适定的，若存在常数0>ε和0>K ，对任意满足条件εδ≤及εη≤∞)(x 的δ和)(x η，常微分方程初值问题⎪⎩⎪⎨⎧+=<<+=δηa a z b x a x z x f dx dz)(),(),(（1.4）存在唯一解)(x z ，且}.{)()(δη+≤-∞∞K x z x y适定问题的解)(x y 连续依赖于（1.1）右端的),(y x f 和初值α。
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（IVP）即为一种最常见的微分方程求解问题，其求解方法有多种，本文将对常微分方程初值问题的数值解法进行较为详细的介绍。
一、欧拉法欧拉法是最基本的一种数值解法，它采用泰勒级数展开并截断低阶项，从而获得一个差分方程近似求解。
具体来讲，设 t 为独立变量，y(t) 为函数 y 关于 t 的函数，方程为：$$y'(t) = f(t, y(t)), \qquad y(t_0) = y_0$$其中 f(t,y(t)) 为已知的函数，y(t_0) 为已知的初值。
将函数 y(t) 进行泰勒级数展开：$$y(t+h) = y(t) + hf(t, y(t)) + O(h^2)$$其中 h 表示步长，O(h^2) 表示其他高阶项。
为了使误差较小，一般取步长 h 尽可能小，于是我们可以用欧拉公式表示数值解：$$y_{n+1} = y_n + hf(t_n, y_n), \qquad y_0 = y(t_0)$$欧拉法的优点是容易理解和实现，但是由于截取低阶项且使用的单步法，所以误差较大，精度较低，在具体应用时需要慎重考虑。
二、龙格-库塔法龙格-库塔法（Runge-Kutta method）是一种多步法，比欧拉法更加精确。
龙格-库塔法的主要思想是使用不同的插值多项式来计算近似解，并且将时间步长分解，每次计算需要多次求解。
以下简要介绍二阶和四阶龙格-库塔法。
二阶龙格-库塔法将时间步长 h 分解成两步 h/2，得到近似解表达式：$$\begin{aligned} k_1 &= hf(t_n, y_n)\\ k_2 &= hf(t_n+h/2,y_n+k_1/2)\\ y_{n+1} &= y_n+k_2+O(h^3)\\ \end{aligned}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。
第八章常微分方程的数值解法
y( xn1 )
15
Euler法的收敛性
称初值问题(8.1.1)的数值解法是收敛的，如:
h0 ( n )
lim yn y ( x)
其中： x xn x0 nh , x [ x0 , b]
16
例考察以下初值问题Euler法的收敛性
dy y dx y (0)＝y0 ( 0)
★
可得： h (k ) ( k 1) y y | f ( xn 1 , yn ) f ( x , y 1 n 1 n 1 ) | 2 hL ( k ) hL k 1 (1) ( k 1) (0) | yn 1 yn 1 | ( ) | yn 1 yn 1 | 2 2 hL k 1 ( k 1) 从而 : lim( ) 0 , 故有 lim yn 1 y n 1 。 k 2 k
★
由y0＝y( x0 ), 假定yn＝y( xn ), 往证：
y0 yn 1 y ( xn 1 ) xn 1； x0
14
证明
yn yn1 yn hf ( xn , yn ) yn h xn 1 1 yn (1 h ) y( xn )(1 h ) xn xn y0 y0 1 xn (1 h ) ( xn h) x0 xn x0 y0 xn 1 x0
8
局部截断误差
假设第n步在点xn的值计算没有误差,即yn y( xn ), 由单步法计算出yn1 , 则
Tn1 y( xn1 ) yn1 称为点xn1上的局部截断误差.
从初值y( x0 ) y0出发，由单步法显式或隐式 逐步计算，得xn 1的值yn 1 , 则
n1 y( xn1 ) yn1
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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 法（i.4）。
如果用右矩形公式，便得到下面的：111Euler (,), 0,1,,1n n n n u u h f t u n N +++=+=-隐式方法： （i.6） 类似地，如果用梯形公式，就得到11111E u l e r [(,)(,)], 0,1,,12n n n n n n hu u f t u f t u n N +++++=++=-改进的方法 (i.7) 当(,)f t u 关于u 是非线性函数的时候，不能由（i.6）和(i.7)从n u 直接算出1n u +，称这一类方法为隐式，通常采用某种迭代法求解。
例如，将一般的隐式方法写成11(,,)n n n n u F t u u ++= (i.8) 则可以利用如下的迭代法由n u 算出1n u +：11101(,,), 0,1, k k n n n n n n u F t u u k u u ++++⎧==⎨=⎩ (i.9)关于k 的迭代通常只需进行很少几步就可以满足精度要求了。
为了避免对隐式方法进行迭代的麻烦，比如说对于改进的Euler 方法（i .7），可以采用某种预估法近似算出11(,)n n f t u ++，然后再用（i .7）作校正，这就导致所谓预估校正法。
下面给出一个例子：1111111(,)2(,)()2n n n n n n n n n n n n n u f t u u u hu u f t u h u u u u +-+++++'⎧=⎧⎪⎪'=+⎨⎪⎪⎪'=⎨⎩⎪⎪''=++⎪⎩预估： 校正： （i.10） 这是一个多步法，即计算节点1n t +上的近似值1n u +时，除了用到前一点的近似值n u 之外，还要用到1n u -，甚至可能用到2,n u -。
而用前面的各种Euler 法计算节点1n t +上的近似值1n u +时，只用到n u ，因此称之为单步法。
下面给出另一个多步法的例子。
在区间2[,]n n t t +上积分（i.1a ），得22()()(,())n n t n n t u t u t f t u t dt ++=+⎰用Simpson 公式(即把被积函数看作二次函数)近似计算积分，便得到212M i l n e (4), 0,1,,23n n n n n hu u f f f n N +++=+++=-法： （i.11） 用多步法（i.10）或（i.11）计算时，必须先用某种单步法由0u 计算出1u ，称为造表头。
然后再逐次算出23,,,N u u u 。
一般说来，多步法比Euler 法等简单的单步法精度要高一些。
下面我们讨论一类所谓Ronge-Kutta 法。
他们是单步法，但是其精度可以与多步法比美。
最常用的是下面的标准Ronge-Kutta 法和Gear 法：112234311234(,)(,)22Lunge-Kutta (,)22(,)(22)6n n n n n n n n n K f t u K h K f t u K h K f t u K f t h u K h u K K K K +=⎧⎪⎪=++⎪⎪⎪=++⎨⎪⎪=++⎪⎪=+++⎪⎩标准法： （i.12）(i.13) 从几何上， Ronge-Kutta 法可以粗略地解释为：在区间1[,]n n t t +中选取若干个点(可以重复)1211n k k n t t ττττ-+=<≤≤≤=，仅仅利用在区间1[,]n n t t +内可以得到的所有信息，依次给出函数(,())f t u t 在这些点上尽可能精确的的近似值12,,,k K K K ，然后把它们组合起来，尽可能精确地近似计算（i.5）中的积分。
Ronge-Kutta 法一般的构造方式如下。
选定常数k ，令1(,)n n K f t u =22211(,)n n K f t h u K αβ=++ 33311322(,)n n K f t h u K K αββ=+++1122,11(,)k n k n k k k k k K f t h u K K K αβββ--=+++++ 11122()n n k k u u h K K K ωωω+=++++ （i.14） 选取这些待定常数,,i ij m αβω的原则是：将（i.14）在(,)n n t u 作Taylor 展开，然后按照h 的幂重新整理，使得231123112!31n n u u h h h γγγ+=++++ （i.15） 与微分方程（i.1a ）的解()u u t =在n t t =处的Taylor 展式23111()()2!3!n n n n n u t u t f h f h f h +'''=++++ （i.16） 有尽可能多的项重合，即要求123, , , n n n f f f γγγ'''=== 这里(,)(,)(,)(,), n n n n n n t u t u df t u f f t u f dt ='==等等。
按照（i.14）构造出的都是显式Runge-Kutta 方法, 在每一个i K 的表达式中只出现n u 。
如果允许在某一个i K 的表达式中出现1n u +，则可以导出隐式Runge-Kutta 方法。
i.2 常微分方程组与高阶常微分方程先来考虑下面的常微分方程组初值问题11122212121100(,,,,)(,,,,)(,,,,)(0),,(0)m m mm m m m du f t u u u dt du f t u u u dt du f t u u u dtu u u u ⎧=⎪⎪⎪=⎪⎪⎨⎪⎪=⎪⎪⎪==⎩ (i.17)利用向量记号,上式可以改写为0()(,)(0)t t u '=⎧⎨=⎩u f u u (i.18) 上节中各方法都可以直接应用到常微分方程组(i.18)。
例如，Euler 方法成为1(,), 0,1,,1n n n n h t n N +=+=-u u f u再来考虑高阶常微分方程111211(,,,,), 0(0),,(0),(0)m m m m m m m d u du d u f t u t dt dt dt d u du v v u v dt dt ----⎧=>⎪⎪⎨⎪===⎪⎩ (i.19) 这时,可以令1121(,,,)(,,)m m m du d u u u u u dt dt --=≡u (i.20) 231212(,,,)(,,,,(,,,,))m m m f f f u u u f t u u u ==f (i.21) 012(,,,)m v v v =u (i.22)于是可以把高阶常微分方程(i.19)化成一阶常微分方程组(i.18)。
i.3 收敛性与稳定性截断误差 粗略地说，截断误差可以定义为将微分方程解带入到差分方程后得到的误差，代表了微分方程与差分方程之间的误差。
例如，由Taylor 展式和微分方程(i.1a)得到221(,())()()()()()(,())2!2!nn n n n n n n t h h df t u t u t u t hu t u u t hf t u t dt ττ+='''=++=++ 其中n τ是区间1[,]n n t t +上某个常数。
与Euler 法1(,)n n n n u u hf t u +=+相比较， 定义余项2(,())2!nt h df t u t dt τ=为Euler 法的截断误差，它关于h 是2阶的，记为2()h O 。
将上节中讨论的单步法写成一般形式1(,,)0n n F u u h +=，则可以定义截断误差为1((),(),)n n F u t u t h +。
对于每一个差分法在适当的点（例如n t t =，1n t +或者11()2n n t t ++）作Taylor 展开，就可以得到截断误差的阶。
对多步法可以类似处理。
上节中各方法的截断误差阶分别为：相容性 一个差分方法称为相容的，如果其截断误差至少是一阶的。
收敛性 实用中我们更关心的是近似解的收敛性，即0h →时，是否有()n n u u t →。
在适当的条件下，例如步长h 足够小，右端函数f 和解u 足够光滑等等，可以证明以上讨论的各方法都是收敛的，并且有估计式()pn n u t u Ch -≤ （i.23） 其中常数C 与u 和f 有关，与h 和n 无关；阶数p 等于截断误差的阶数减去1。