常微分方程初值问题的数值解法
常微分方程初值问题数值解法
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
微分方程初值问题的数值解法
积分法:
yk 1 yk h f ( xk , yk ) y ( x0 ) y0
积分项利用矩形公式计算
(1) y( xk 1 ) y( xk )
xk 1
xk
f (t , y(t ))dt
(★)
xk 1
xk
f (t , y(t ))dt h f ( xk , yk ) y( xk 1 ) y( xk ) h f ( xk , yk )
引言
初值问题的数值解法:求初值问题的解在一系列节点的值 y ( xn )的近似值 yn 的方法.本章数值解法的特点:都是采用“步进 式”,即求解过程顺着节点排列的次序一步步向前推进. 常微分方程初值问题: dy f ( x, y ), x [a, b] dx y ( x0 ) y0
替 f (x , y)关于 y 满足Lipschitz条件. 除了要保证(1)有唯一解外,还需保证微分方程本身是稳定的,即 (1)的解连续依赖于初始值和函数 f (x , y). 也就是说, 当初始值 y0 及函数 f (x , y)有微小变化时, 只能引起解的微小变化.
注: 如无特别说明,总假设(1)的解存在唯一且足够光滑. 在 f 连续有界, 则 f (x , y)对变量 y 可微的情形下, 若偏导数 y 可取L为
也称折线法 x
2. 梯形法
若采用梯形公式计算(★)中的积分项,则有 h y ( xk 1 ) y ( xk ) [ f ( xk , y ( xk )) f ( xk 1 , y ( xk 1 ))] 2 h yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk 1 )] 2 称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做 法: (0) (0) 先用Euler法求出初值 yk ,1 即 ,将其代入梯形公式 yk 1 yk h f ( xk , yk ) 的右端,使之转化为显式公式,即 h ( l 1) (l ) yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk (☆ ) 1 )] 2
第五章 常微分方程初值问题数值解法
则有
yn 1 yn hf ( xn , yn )
( 5.2 ) Euler格式
例5.1 用Euler格式解初值问题
2x y y y y (0) 1
取步长h=0.1.
(0 x 1)
Euler格式的具体形式为
y n 1 y n hf ( x n , y n ) 2 xn yn 0.1( yn ) yn 0.2 xn 1.1 yn yn
计算公式的精度 常以Taylor展开为工具来分析计算公式的精度. 为简化分析，假定yn是准确的，即在 yn y ( xn ) 的前提下估计误差 y ( xn 1 ) yn 1 Euler格式的局部截断误差 由 从而 局部截断误差
f ( xn , yn ) f ( xn , y ( xn )) y '( xn ) y ( xn 1 ) yn 1 y ( xn 1 ) ( yn hf ( xn , yn )) y ( xn 1 ) y ( xn ) hy '( xn )
y ( xn ), y ( xn 1 ), 的近似值 y1 , y2 , , yn , yn 1 ,
相邻两个节点的间距 h xi 1 xi 称为步长，步 长可以相等，也可以不等.本章总是假定h为定数， 称为定步长，这时节点可表示为
xn x0 nh , n 0,1, 2,
由f ( xn 1 , yn 1 ) f ( xn 1 , y ( xn 1 )) f y ( xn 1 , )( yn 1 y ( xn 1 )) f ( xn 1 , y ( xn 1 )) y '( xn 1 )（在xn点Taylor展开） h2 y '( xn ) hy ''( xn ) y '''( xn ) ... 2 3 2 h h 因此yn 1 y ( xn ) hy '( xn ) y ''( xn ) y '''( xn ) 2 4 h f y ( xn 1 , )( yn 1 y ( xn 1 )) 2 h2 h3 y ( xn 1 ) y ( xn ) hy '( xn ) y ''( xn ) y '''( xn ) 2 3!
第9章 常微分方程初值问题数值解法
数值分析
第9章 常微分方程初值问题数值解法
《常微分方程》中介绍的微分方程主要有:
(1)变量可分离的方程 (2)一阶线性微分方程(贝努利方程) (3)可降阶的一类高阶方程 (4)二阶常系数齐次微分方程 (5)二阶常系数非齐次微分方程 (6)全微分方程 本章主要介绍一阶常微分方程初值问题的数值解法。
进一步: 令
y n1 y n
xn 1 xn
y n 1 y( x n 1 ) , y n y( x n )
f ( x , y( x ))dx h f ( x n , y n )
宽
9
高
实际上是矩形法
数值分析
第9章 常微分方程初值问题数值解法
(3)
用Taylor多项式近似并可估计误差
解决方法：有的可化为显格式，但有的不行 18
数值分析
第9章 常微分方程初值问题数值解法
与Euler法结合，形成迭代算法 ，对n 0，2， 1，
( yn0 )1 yn hf x n , yn ( k 1) h ( yn1 yn f x n , yn f x n1 , ynk )1 2
7
数值分析
第9章 常微分方程初值问题数值解法
建立数值解法的常用方法
建立微分方程数值解法,首先要将微分方程离散 化. 一般采用以下几种方法: (1) 用差商近似导数
dy yx yx x x dx x y
n 1 n n 1 n
n
,
n
进一步: 令
yn1 y( xn1 ) , yn y( xn )
由 x0 , y0 出发取解曲线 y y x 的切线（存在！），则斜率
第9章 常微分方程初值问题数值解法
oa
b
a f ( x)dx (b a) f (b)
中矩形公式
b
ab
a f ( x)dx (b a) f ( 2 )
计算方法
梯形公式
bx
右矩形公式 中矩形公式 左矩形公式
§ 欧拉方法几何意义
y y y(x)
y0 y1 y2 0 x0 x1 x2
计算方法
x
§ 隐式欧拉方法
➢隐式欧拉法 /* implicit Euler method */
初 值 问 题 的 解 必 存 在 且唯 一 。
计算方法
§9.1 引言
三. 数值解法含义
所谓数值解法, 就是设法将常微分方程离散化, 建 立差分方程, 给出解在一些离散点上的近似值。
微分方程的数值解: 设方程问题的解y(x)的存在区 间是[a,b], 令a= x0< x1<…< xn =b, 其中hk=xk+1-xk, 如是等距节点h=(b-a)/n, h称为步长。
yi1 yi1 2h f ( xi , yi ) i 1, ... , n 1
计算方法
预估-校正法
三. 预估 — 校正法
/* predictor-corrector method */
方法 显式欧拉 隐式欧拉 梯形公式
中点公式
简单
稳定性最好
精度提高
精度低
精度低, 计算量大
计算量大
精度提高, 显式
在x0 x X上的数值解法。
四. 误差估计、收敛性
和稳定性
计算方法
§9.2 简单的数值方法与基本概念
一. 欧拉（Euler）格式
设 节 点 为xi a ih (i 0,1,2 , n) 方 法 一 ：Taylor展 开 法
第6章常微分方程初值问题的解法
ykh 2 k[ (ykx k 1 ) ( yk 1x k 1 1 )]
yk11 29 1yk1k05110
预估－校正Euler方法：
y k 1 0 .90 y k 5 0 .00 k 9 0 .1 5
20
Euler方法
xk
yk
yk y(xk)
0.0 1.000000
0.0
梯形方法
yk
yk y(xk)
1.000000
0.0
续
预估－校正方法
yk
yk y(xk)
1.000000
0.0
0.1 1.000000 0.2 1.010000
4.8×10－3 8.7×10－3
1.004762 1.018594
y(0) 1
其解析解为： y1xe-t2dt x[0,1] 0 很难得到其解析解
4
例如：
y=x+y , x[0,1]
y(0) 1
其解析解为 yx12ex
只有一些特殊类型的微分方程问题能够得到用解析表达式 表示的函数解，而大量的微分方程问题很难得到其解析解。
因此，只能依赖于数值方法去获得微分方程的数值解。
例如：
y=x+y , x[0,1]
y(0) 1
其解析解为：yx12ex
3
但是, 只有一些特殊类型的微分方程问题能够得到用解析 表达式表示的函数解，而大量的微分方程问题很难得到其解 析解。
因此，只能依赖于数值方法去获得微分方程的数值解。
例如：
y =e-x2 ,
x[0,1]
7.5×10－5 1.4×10－4
常微分方程初值问题的数值解法
第七章 常微分方程初值问题的数值解法--------学习小结一、本章学习体会通过本章的学习，我了解了常微分方程初值问题的计算方法，对于解决那些很难求解出解析表达式的，甚至有解析表达式但是解不出具体的值的常微分方程非常有用。
在这一章里求解常微分方程的基本思想是将初值问题进行离散化，然后进行迭代求解。
在这里将初值问题离散化的方法有三种，分别是差商代替导数的方法、Taylor 级数法和数值积分法。
常微分方程初值问题的数值解法的分类有显示方法和隐式方法，或者可以分为单步法和多步法。
在这里单步法是指计算第n+1个y 的值时,只用到前一步的值，而多步法则是指计算第n+1个y 的值时，用到了前几步的值。
通过对本章的学习，已经能熟练掌握如何用Taylor 级数法去求解单步法中各方法的公式和截断误差，但是对线性多步法的求解理解不怎么透切，特别是计算过程较复杂的推理。
在本章的学习过程中还遇到不少问题，比如本章知识点多，公式多，在做题时容易混淆，其次对几种R-K 公式的理解不够透彻，处理一个实际问题时，不知道选取哪一种公式，通过课本里面几种方法的计算比较得知其误差并不一样，，这个还需要自己在往后的实际应用中多多实践留意并总结。
二、本章知识梳理常微分方程初值问题的数值解法一般概念步长h ，取节点0,(0,1,...,)n t t nh n M =+=，且M t T ≤，则初值问题000'(,),()y f t y t t Ty t y =≤≤⎧⎨=⎩的数值解法的一般形式是1(,,,...,,)0,(0,1,...,)n n n n k F t y y y h n M k ++==-@显示单步法7.2.1 显示单步法的一般形式1(,,),(0,1,...,1)n n n n y y h t y h n M ϕ+=+=-定理7.2.1 设增量函数(,,)n n t y h ϕ在区域00{(,,)|,||,0}D t y h t t T y h h =≤≤<∞≤≤内对变量y 满足Lipschitz 条件，即存在常数K ，使对D 内任何两点1(,,)t u h 和2(,,)t u h ，不等式1212|(,,)(,,)|||t u h t u h K u u ϕϕ-≤-成立，那么，若单步法的局部截断误差1n R +与1(1)p h p +≥同阶，即11()p n R O h ++=，则单步法的整体截断误差1n ε+与p h 同阶，即1()p n O h ε+=。
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解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乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
常微分方程初值问题的数值积分法
y( xn1) y( xn ) hfh ( xn , y( xn )) Rh ( xn ),
并且当 h 0时，
yh,0
y(x0 )
h1 max xIh
Rh (x)
0,
(7.4.2)
则称(7.4.1) 式为初值问题(7.1.1)的一个相容近似 ，
或称此格式满足相容条件即(7.4.2)式。
若
Rh (xn ) O(h p1), yh,0 y(x0) O(h p )(h 0)
a x0 x1 L xN1 xN b,
令 hn xn1 xn，称为积分网格的步长。
用y0 , y1, y2 , , yN 表示精确解 y(x)在节点 x0 , x1, x2 , , xN 上函数值 y( x0 ), y( x1), y( x2 ), , y( xN )的近似值。
对给定的数值积分法，各个 yn 是按某一递推算法确 定的。若在计算 yn1 时只用到已求出的 y0, y1,L , yn中的 yn ，而无须使用其余值 y0, y1,L , yn1 中的任何一个， 则称此法为单步法，否则，称之为多步法。
xn
的右端积分中用梯形公式，则得
yn1
yn
h 2
[
f
(
xn
,
yn
)
f (xn1, yn1)],
n 0,1,L
称该递推公式为梯形方法。
梯形公式
b f (x)dx (b a) ( f (a) f (b))
a
2
梯形方法
yn1
yn
h 2
[
f
( xn
,
yn
)
f (xn1, yn1)],
n 0,1,L
7.2 几个简单的数值积分法
第7章-常微分方程初值问题的数值解法
由 点 斜 式 写 出 切 线 方 程 ：
dy yy0(xx0)dx(x0,y0) y0(xx0)f(x0,y0)
等 步 长 为 h ， 则 x x 0 h ， 可 由 切 线 算 出 y 1
y1y0hf(x0,y0)
2021/4/9
10
其 中 ， y n 1 是 当 y n y (x n )(精 确 解 )时 由 E u le r法 求 出 的 值 ， 即 y n 无 误 差 。
将 y (x n 1 )在 x n 点 T a y lo r展 开 ：
y (x n 1 ) y (x n h ) y (x n ) h f(x n ,y (x n )) (h 2 x2 n y ()xn1)
ax0 x1 x2 xn b 在节点上用离散化方法将连续型微分方程 转化成离散型代数方程即差分方程来求解。
具 体 作 法 ： 利 用 y(x0)求 出 y(x1)的 近 似 值 y1， 再 由
y1求 出 y2，， 直 到 求 出 yn为 止 。 该 算 法 称 为 步 进
式 或 递 推 式 算 法 。
积 分 用 梯 形 公 式 ， 且 令 ： yn1y(xn1),yny(xn) 则 得 ： yn1ynh 2(f(xn,yn)f(xn1,yn1))
R n 1 y (x n 1 ) y n 1 1 h 2 3y () x n x n 1
与Euler法结合，形成迭代算法，对n0， 1， 2，
yn(0)1ynhf(xn,yn) (75)
y 2021/4/9
(k1) n1
ynh2(f(xn,yn)f(xn1,yn(k1))k0,1,2,
12
2. Euler方法的截断误差
计算方法课件第八章常微分方程初值问题的数值解法
整体截断误差与局部截断误差的关系
定理：如果f(x，y)满足李普希兹（Lipschitz）条件
f(x ,y 1 )f(x ,y 2) L y 1y 2
且局部截断误差有界：
|R n|1 2h2M 2
(n1,2, )
则Euler法的整体截断误差n满足估计式：
ne(ba)L 0h 2L M 2(e(ba)L1)
分光滑。初值问题的解析解(理论解）用 y(x表n ) 示, 数值解法的精确解用 y表n 示。
常微分方程数值解法一般分为：
（1）一步法：在计算y n 1 时，只用到x n 1 ,x n和 y，n 即前一步的值。
（2）多步法：计算 y n 1 时，除用到 x n 1 ,x n 和 y n 以外，还要用 x n p 和 y n p (p1 ,2 k;k0) ,即前
其中L为李普希兹常数，b-a为求解区间长度，
M2 mayx(x) 。 axb
证明参见教材。
Remark：该定理表明，整体截断误差比局部截 断误差低一阶。对其它方法，也有类似的结论。
收敛性与稳定性
收敛性定义：如果某一数值方法对于任意固定的
xn=x0+nh，当h0(同时n )时有yn y(xn)，
则称该方法收敛。 稳定性定义 定义 用一个数值方法，求解微分方程初值问 题时，对给定的步长h>0，若在计算 y n 时引入 误差 （n 也称扰动），但由此引起计算后面的 ynk(k1,2, )时的误差按绝对值均不增加，则 称这个数值方法是稳定的。
一般的显式rk方法可以写成型钢截面只需少量加工即可用作构件省工省时成本低但型钢截面受型钢种类及型钢号限制难于完全与受力所需的面积相对应用料较多其中为常数选取这些常数的原则是要求第一式的右端在处泰勒展开后按h型钢截面只需少量加工即可用作构件省工省时成本低但型钢截面受型钢种类及型钢号限制难于完全与受力所需的面积相对应用料较多上述公式叫做n级的rungekutta方法其局部截断误差为显然euler法是一级一阶rk方法
第一章常微分方程初值问题数值解法
（1.2.3）
其中rn,k(t)为插值余项。 代到（1.2.2）式中得
u ( tn +1 ) = u ( tn ) +
舍去余项 并用uj代替u(tj)即得
∫
tn+1 tn
Ln , k ( t ) dt + ∫ t rn , k ( t ) dt
tn+1
n
(1.2.4 (1.2.5)
Rn , k = ∫
⎡ ∑ ⎣α u
j =0
j =0
αk ≠ 0
（1.2.1）
j n+ j
⎤ − hβ j f n + j ⎦ = 0（数值解满足的差分方程）
因此称（1.2.1）为多步法 或 k-步法。 又因为(1.2.1)关于 u n + j , f n + j 是线性的，所以称为线性多步法。 为使多步法的计算能够进行，除给定的初值u0 外，还要 知道附加初值u1,u2,…,uk-1 ，这可用其它方法计算。 若 β k = 0 则称(1.2.1)是显式的； 若 β k ≠ 0 则方法（1.2.1）是隐式的。 例如，一般线性二步法可写成：
f ( t , u ( t ) ) = Ln , k +1 ( t ) + rn , k +1 ( t )
其中rn,k+1(t)为插值余项。 同理即
un +1 = un + h ∑ bk +1i f ( tn −i +1 , un −i +1 )
i =0
k +1
其中
bk +1i
=∫ ∏
−1
j =0 j ≠i
0
k +1
9、常微分方程初值问题数值解法
( k +1) yn +1
xn +1 ∫xn
− yn +1 |≤
hL 2
|
(k ) yn +1
− yn +1 |,
( 只要 hL < 1，则( 2.8)的ynk +1)收敛到(2.7)的yn +1. +1 2
三、单步法的局部截断误差与阶
一阶常微分方程初值问题(1.1)(1.2)的单步法的一般形式 yn +1 = yn + hϕ ( xn , yn , yn +1, h).
clear x=0,yn=1 %初始化 for n=1:10 yp=yn+0.1*(yn-2*x/yn); %预测 x=x+0.1; yc=yn+0.1*(yp-2*x/yp) ; yn=(yp+yc)/2 %校正 end
( 2.2)
作业： 作业：P381, 1, 2(1).
龙格—库塔 库塔(Runge-Kutta)法 §3 龙格 库塔 法
进一步 y ( xn +1 ) = y ( xn ) + ∫
xn +1 xn
f ( x, y ( x))dx,
(Байду номын сангаас.3)
∫
⇒ 其中
xn +1 xn
f ( x, y ( x))dx ≈ h ∑ ci f ( xn + λi h, y ( xn + λi h)).
yn +1 = yn + hϕ ( xn , yn , h),
i =1
r
(3.4) (3.5) 欧拉法r = 1, p = 1.改进
常微分方程初值问题的的数值解法
本章讨论常微分方程初值问题的数值解法
2
考虑一阶常微分方程的初值问题
⎧ dy ⎪ = f ( x, y ) ⎨ dx ⎪ ⎩ y (a ) = y0
x ∈ [a, b]
只要 f (x, y) 在[a, b] × R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使对任意x∈[a, b] ，和y1， y2 ∈ R1 都有 | f ( x, y1) − f ( x, y2 ) | ≤ L| y1 − y2 | 在唯一解。 成立， 则上述问题存
⎧ ⎪ ⎨ ⎪ ⎩ y n +1 = yn + hf ( xn , yn ), h yn +1 = yn + [ f ( xn , yn ) + f ( xn +1 , y n +1 )] 2
改进的Euler方法：y0=1,
y1=y0+hf (x0, y0) =1.1, y1=1+01./2 ×[(1−2 ×0/1)+(1.1−2 ×0.1/1.1)] =1.095909, …… y11=…… y11=1.737869.
1 yn +1 = yn + h[ f ( xn , yn ) + f ( xn +1 , yn +1 )] 2
12
称之为梯形公式。这是一个隐式的计算公式，欲求的yn+1需 解一个方程。
3．截断误差
定义 在假设 yn = y(xn)，即第 n 步计算是精确的前提下，考 虑的截断误差 εn+1 = y(xn+1) − yn+1 称为局部截断误差
⎧ y n +1 = y n + k1 ⎨ ⎩k1 = hf ( xn ,y n )
第六章常微分方程的数值解法
第六章常微分方程的数值解法第六章常微分方程的数值解法在自然科学研究和工程技术领域中，常常会遇到常微分方程的求解问题。
传统的数学分析方法仅能给出一些简单的、常系数的、经典的线性方程的解析表达式，不能处理复杂的、变系数的、非线性方程，对于这些方面的问题，只能求诸于近似解法和数值解法。
而且在许多实际问题中，确确实实并不总是需要精确的解析解，往往只需获得近似的解或者解在若干个点上的数值即可。
在高等数学课程中介绍过的级数解法和逐步逼近法，能够给出解的近似表达式，这一类方法称为近似解法。
还有一类方法是通过计算机来求解微分方程的数值解，给出解在一些离散点上的近似值，这一类方法称作为数值方法。
本章主要介绍常微分方程初值问题的数值解法，包括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]上是连续可微的。
常微分方程初值问题数值解法
y ( xn ) 1.0954 1.1832 1.2649 1.3416
1.4351 1.4142
初值问题（2.2）有解 y 1 2 x ，按这个解析式子 算出的准确值 y ( xn ) 同近似值 y n一起列在表9-1中，两者 相比较可以看出欧拉方法的精度很差. 还可以通过几何直观来考察欧拉方法的精度. 假设 yn y ( xn ) ，即顶点 Pn 落在积分曲线 y y ( x) 上，那么，按欧拉方法做出的折线 Pn Pn 1 便是 y y ( x) 过点 Pn 的切线（图9-2）.
yn 1 yn hf ( xn 1 , yn 1 ),
（2.5）
称为后退的欧拉法. 后退的欧拉公式与欧拉公式有着本质的区别，后者是 关于 yn 1 的一个直接的计算公式，这类公式称作是显式的；
即
yn1 yn hf ( xn , yn ). y0 r）公式. 式（2.1）可逐步算出
y1 y0 hf ( x0 , y0 ), y2 y1 hf ( x1 , y1 ),
例1
求解初值问题
2x y y y y (0) 1. (0 x 1),
n ( xn , xn 1 ).
在 yn y ( xn ) 的前提下， f ( xn , yn ) f ( xn , y( xn )) y( xn ) . 于是可得欧拉法（2.1）的公式误差
y ( xn 1 ) yn 1 h2 h2 y( n ) y( xn ), 2 2
（1.3）
y y ( x) 理论上就可以保证初值问题（1.1），（1.2）的解 存在并且唯一.
1
所谓数值解法，就是寻求解 y ( x) 在一系列离散节点
计算方法 常微分方程初值问题数值解法 Euler公式 龙格 库塔法
第12次 常微分方程初值问题数值解法
内容
1. 常微分方程初值问题解的存在性定理 2. Euler公式 3. 梯形公式 4. 两步Euler公式 5. 欧拉法的局部截断误差 6. 改进型Euler公式 7. 龙格-库塔法 8. 算法实现
常微分方程初值问题 解的存在性定理
欧拉法的局部截断误差
9.2.4． 欧拉法的局部截断误差
衡量求解公式好坏的一个主要标准是求解公式的精 度, 因此引入局部截断误差和阶数的概念。
定义9.1 在yi准确的前提下, 即 yi y(xi)时, 用数值
方法计算yi+1的误差：
R i yi( 1 ) x y i 1
称为该数值方法计算时yi+1的局部截断误差。
通常表示成下列平均化形式
yi1的近似 y p y i hf(x i , y i )
yi1的近似 y c y i hf(x i 1 , y p )
yi1的更好 的近似
y i1
1 2
(y p
yc)
i 0, 1, 2 … , n 1
(9.12)
例9.2 用改进欧拉法解初值问题
y
y
2x y
f( y 1 ) x f,( y 2 ) x L y , 1 y 2 , y 1 ,y 2 R
则方程( 9.1 ) 在a, b上存在唯一的连续可微分的 解的解 y=y(x) 。
推论：如果函数f(x,y)对y的偏导数
f(x, y
y)
在带形
区域 R { a x b ,y - }
内有界。
两步欧拉公式的局部截断误差为：
y(i1 x )yi1O3 ()h
从而两步欧拉公式的阶数是2.推导过程省略。
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（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}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第九章 常微分方程初值问题的数值解法第一部分 内容提要一、数值解的一般概念常微分方程初值问题00'()(,)()y x f x y y x y =⎧⎨=⎩的数值解是指通过一定的近似方法得出准确解()y y x =在一列离散点012,,,,,n x x x x 上的近似值012,,,,,n y y y y 。
数值解的特征是步进式，即()y x 在1n x +点的近似值1n y +是由1,,n n x x -等若干点处的近似值1,,n n y y -的信息给出的递推公式。
若1n y +依赖于前面k 步的值11,,,n n n k y y y --+，则称为k 步法；1k =称为单步法。
利用()y x 在11,,,n n n k x x x --+的精确解11(),(),,()n n n k y x y x y x --+借助某种算法计算出1n y +，则称11()n n y x y ++-为该方法的局部截断误差。
如果一个算法的局部截断误差是1()p O h +，则称该方法是p 阶的；而利用数值解11,,,n n n k y y y --+得到的1n y +与微分方程的精确解之差11()n n y x y ++-称为整体截断误差，即是该数值方法的误差。
对于固定的0x x >，取0x x h n-=，用某种算法得到n y ，如有lim ()n h y x y →-=0，则称该方法是收敛的。
注意，因x 是固定的，随着0h →，数值解的步数n →∞。
在实际计算时由于舍入误差不可避免，实际得到数值解是n y ，稳定性即研究n n y y -是否随着计算步骤n 的增加而增加。
通常所提的稳定性是通过模型方程(0)y y λλ'=<来讨论的。
若当某一步n y 有舍入误差时，在以后的计算中误差不会逐步扩大，则称这种稳定性为绝对稳定性。
二、简单单步法及其收敛性、稳定性Euler 法1(,)n n n n y y hf x y +=+的局部截断误差为2()O h ，整体截断误差为()O h ，即一阶收敛。
对于模型问题(0)y y λλ'=<，当20h λ<<-时，Euler法是数值稳定的。
隐式Euler 法111(,)n n n n y y hf x y +++=+的误差与Euler 法相同，但是无条件稳定：即对任意步长0h >，隐式Euler 法都是稳定的。
梯形法[]111(,)(,)2n n n n n n hy y f x y f x y +++=++的误差比Euler 法高一阶，也是无条件稳定的。
改进Euler 法[]11(,)(,(,))2n n n n n n n n hy y f x y f x y hf x y ++=+++是一种预测－校正方法：1(,)p n n n n y y hf x y +=+——Euler 法预测111(,)(,)2c p n n n n n n h y y f x y f x y +++⎡⎤=++⎣⎦——梯形法校正 它保持了梯形法的误差阶数，但不是无条件稳定的。
三、龙格-库塔方法龙格-库塔类算法采用区间[]1,n n x x +内若干点的斜率的加权平均来近似整个区间的平均斜率，一般形式为11111(,)(,),2,,s n n i ii n n i i n i n i ij j j y y h K K f x y K f x c h y c h a K i sλ+=-=⎧=+⎪⎪⎪=⎨⎪⎪=++=⎪⎩∑∑如经典的4级(4)s =4阶（局部截断误差为5()O h ）Runge-Kutta 公式为112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （1）四、线性多步法线性k 步法具有一般形式1011110111()n n n k n k n n k n k y y y y h f f f αααβββ+---++-+=+++++++ （2）00≠β为隐式公式；00=β为显式公式。
构造多步法公式有基于数值积分和Taylor 展开两种途径。
多步法（2）的局部截断误差为111110()()(,())k kn n j n j j n j n j j j T y x y x h f x y x αβ-++--+-+===--∑∑利用原微分方程后，成为11110()()().k kn n j n j j n j j j T y x y x h y x αβ-++--+=='=--∑∑因此利用Taylor 公式，分别对()y x 和()y x '作Taylor 展开，可确定线性多步（2）中的待定参数{}10-=k j j α和{}kj j 0=β，使它达到最高阶精度（或指定精度）。
预测－校正格式：不论单步法还是多步法，隐式公式比显式公式的稳定性好，但隐式公式的计算比较困难。
预测－校正格式是用显式公式进行预测，再用隐式公式进行校正。
五、高阶方程和一阶常微分方程组所有应用于一阶常微分方程初值问题的数值方法都可以直接推广到一阶常 微分方程组，只需把公式中的未知函数改为向量形式。
高阶方程总可以降为一阶方程组，进而用数值方法求解。
第二部分 例题精选一、对初值问题y y -='，1)0(=y ，用梯形公式求解，求数值解n y 的表达式（写成步长h 的函数），并证明数值解的收敛性。
分析：微分方程用数值方法离散后即变成差分方程，单步法导出的差分方程通常是关于数值解序列{}n y 的一个递推公式，因此问题变为已知数列首项和递推公式求数列通项公式，收敛性即是数列极限。
解：梯形公式)(211++++=n n n n f f hy y 对本问题n n y f -=于是梯形公式的解1122++--=n n n n y h y h y y 即n n y h h y ⎪⎭⎫ ⎝⎛+-=+221由于10=y ，易得 、、、321,22=⎪⎭⎫⎝⎛+-=n h h y nn显然初值问题的准确解为x e x y -=)(对于固定的点nh x x n ==，准确解在该点的值为n x n e x y -=)(而数值解为hx nn nh h h h y ⎪⎭⎫ ⎝⎛+-=⎪⎭⎫ ⎝⎛+-=22122令n x 固定，上式对0h →取极限得：)(2212222n x h x h h n x y e h h y n n =→⎪⎭⎫ ⎝⎛+-=-⎪⎭⎫⎝⎛+-⋅+-这证明了数值解的收敛性。
二、用隐式单步法（该方法也属于隐式Runge-Kutta 方法）⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++=++=+),()3,3()3(41211211hk y h x f k k h y h x f k k k h y y n n n n n n求解微分方程初值问题y y 5'-=，1)0(=y 时，试推出其绝对稳定区间。
分析：将格式应用于所给方程可导出误差传播方程，从而求出绝对稳定区间 解：记舍入误差为n ε。
该隐式单步法应用于方程y y 5'-=时h y k k h y k n n 531535111+-=⇒⎪⎭⎫ ⎝⎛+-=()n n n n y hh h hy y hk y k 53155053755512+-=++-=+-= n n n n y h h h y h h h y y 1066202553155045421++-=⎪⎭⎫ ⎝⎛+-+-+=+ 从而误差传播方程为n n hh h εε1066202521++-=+解不等式1106620252≤++-hh h 得560≤≤h故，该方法应用到所给方程的绝对稳定区间为]56,0(。
三、用Taylor 展开原理构造形如11011()()n n n n n y y y h f f αββ+--=+++的两步法，试确定系数01,,,αββ使方法具有二阶精度，并推导其局部截断误差主项。
分析：本题考察构造多步法的方法；二阶精度即局部截断误差为3()O h 。
解：由多步法局部截断误差的定义[]111011()()()()().(3)n n n n n n T y x y x y x h y x y x ααββ++--''=---+ 将1()n y x +，1()n y x -和1()n y x -'分别在点n x 展开得：2341()()()()()()()26n n n n n n h h y x y x h y x hy x y x y x O h +''''''=+=++++2341()()()()()()()26n n n n n n h h y x y x h y x hy x y x y x O h -''''''=-=-+-+231()()()()()()2n n n n n h y x y x h y x hy x y x O h -''''''''=-=-++将此三式代入局部截断误差(3)式，我们有23410111(12)()()(1)()(12)()(13)()26n n n n n h h T y x hy x y x y x O h ααββαβαβ+''''''=-++--+-+++-+要使方法具有二阶精度，必需011120,10,120.ααββαβ-=+--=-+=解得01171,,244αββ===-。
此时13138αβ+-=，所以局部截断误差主项为33()8n h y x '''。
四、取0.1h =，试用Euler 法求解初值问题2420,00.2(0)1,(0)0y xyy y x y y '''⎧++=≤≤⎨'==⎩。
分析：先将二阶方程写为一阶方程组，再用Euler 法求解。
解：引进,u y v y '==，并记向量函数u v ⎛⎫= ⎪⎝⎭y ，则原二阶方程变为一阶方程组2(,)42u v x v xuv u '⎛⎫⎛⎫'=== ⎪ ⎪--⎝⎭⎝⎭y f y ，它满足初值条件0(0)1(0)(0)0u v ⎛⎫⎛⎫=== ⎪ ⎪⎝⎭⎝⎭y y 。
将Euler 法1(,)n n n n h x +=+y y f y 应用到该方程组得：1000101(,)0.1020.2h x ⎛⎫⎛⎫⎛⎫=+=+= ⎪ ⎪ ⎪--⎝⎭⎝⎭⎝⎭y y f y211110.20.98(,)0.10.20.4(0.2)20.392h x -⎛⎫⎛⎫⎛⎫=+=+= ⎪ ⎪ ⎪-----⎝⎭⎝⎭⎝⎭y y f y即，用Euler 法求得原二阶方程的数值解为121,0.98y y ==。