常微分方程初值问题数值解法
常微分方程初值问题数值解法
常微分方程初值问题的数值解法在自然科学、工程技术、经济和医学等领域中，常常会遇到一阶常微分方程初值问题：(,),,(),y f x y a x b y a y '=≤≤⎧⎨=⎩ (1) 此处f 为,x y 的已知函数，0y 是给定的初始值。
本章讨论该问题的数值解法，要求f 在区域{(,)|,}G x y a x b y =≤≤<∞内连续，并对y 满足Lipschitz 条件，从而初值问题(1)有唯一的连续可微解()y y x =，且它是适定的。
1 几个简单的数值积分法1.1 Euler 方法（1）向前Euler 公式（显式Euler 公式）10(,),0,1,2,,(),n n n n y y hf x y n y y a +=+=⎧⎨=⎩(2) 其中h 为步长。
由此便可由初值0y 逐步算出一阶常微分方程初值问题(1)的解()y y x =在节点12,,x x 处的近似值12,,y y 。
该公式的局部截断误差为2()O h ，是一阶方法。
（2）向后Euler 公式（隐式Euler 公式）1110(,),0,1,2,,(),n n n n y y hf x y n y y a +++=+=⎧⎨=⎩(3) 这是一个隐格式，也是一阶方法。
这类隐格式的计算比显格式困难，一般采用迭代法求解。
首先用向前Euler 公式提供迭代初值，然后迭代计算：(0)1(1)()111(,),(,),0,1,2,n n n n k k n n n n y y hf x y y y hf x y k +++++⎧=+⎨=+=⎩ (4)1.2 梯形方法1110[(,)(,)],2(),(0,1,2,)n n n n n n h y y f x y f x y y y a n +++⎧=++⎪⎨⎪=⎩= (5) 这也是一个隐格式，是二阶方法。
一般也采用迭代法求解。
迭代公式如下：(0)1(1)()111(,),[(,)(,)],0,1,2,2n n n n k k n n n n n n y y hf x y h y y f x y f x y k +++++⎧=+⎪⎨=++=⎪⎩ (6)1.3 改进的Euler 方法11110(,),[(,)(,)],0,1,2,,2(),n n n n n n n n n n y y hf x y h y y f x y f x y n y y a ++++⎧=+⎪⎪=++=⎨⎪=⎪⎩(7) 为了便于上机编程计算，(7)可改写为110(,),(,),0,1,2,,1(),2(),p n n n cn n p n p c y y hf x y y y hf x y n y y y y y a ++=+⎧⎪=+⎪⎪=⎨=+⎪⎪=⎪⎩(8) 该格式是显式，也是二阶方法。
计算方法 常微分方程初值问题数值解法-Euler公式-龙格-库塔法
[xi , xi 1 ]上积分得，
y(xi 1 ) y(xi )
xi 1
xi
f[x, y(x)]dx
(9.4 )
改用梯形方法计算其积分项，即
xi 1
x i 1 x i [f(x i , y(x i )) f(x i 1 , y(x i 1 ))] 2
xi
f[x, y(x)]dx
0 1 n1 n
… , y(xn ) (未知) 处的函数值 y(x 0 ), y(x1 ),
, yn 的近似值 y 0 , y1 ,…
y=y(x)
a=x0 x1
x2
x3
xn=b
• 相邻两个节点的间距 h xi 1 xi 称为步长，
步长可以相等，也可以不等。
• 本章总是假定h为定数，称为定步长，这时节 点可表示为
第9章 常微分方程初值问题数值解法
§9.1 引言
包含自变量、未知函数及未知函数的导数的方程称 为微分方程。
自变量个数只有一个的微分方程称为常微分方 程。
微分方程中出现的未知函数最高阶导数的阶数 称为微分方程的阶数。 如果未知函数y及其各阶导数
y, y, … , y
(n)
都是一次的,则称其为线性的,否则称为非线性的。
• 如下是一些典型方程求解析解的基本方法 可分离变量法、 常系数齐次线性方程的解法、 常系数非齐次线性方程的解法等。
• 但能求解的常微分方程仍然是很少的，大多数
的常微分方程是不可能给出解析解。例如，一
阶微分方程
y x y
2
2
的解就不能用初等函数及其积分来表达。
• 从实际问题当中归纳出来的微分方程，通常主 要依靠数值解法来解决。 • 本章主要讨论一阶常微分方程初值问题
浅谈常微分方程初值问题数值解法
浅谈常微分方程初值问题数值解法在自然科学、工程技术、甚至社会科学的一些领域中，常常会遇见一阶常微分方程的求解问题：()上述问题，寻求解的具体表达式十分困难，仅对一些特殊形式的才有可能找到解的解析表达式，在大多情况下，初值问题的解不能用初等函数表示出来即使可写出解的解析表达式，但因为这些表达式过于复杂，要计算它在某些点上的函数值也异常困难。
在实际问题中，经常需要的恰是解在某些点上的函数值，因此研究初值问题的数值解法十分必要。
1 常微分方程初值问题的数值解法常微分方程的近似解法大体可分成三大类：一类是图解法和器械法；第二类是解的近似法；第三类是数值解法，即通过离散化的方法直接求出函数在某些点上的近似值，此数值解仅为精确解的近似解。
其基本原理为：一阶常微分方程的初值问题的解是上变量的连续函数，因此求上述问题的数值解，就是在区间上的若干离散点上用离散化的方法将初值问题化成离散变量的相应问题，从而相应问题的解可作为初值问题理论解的近似值。
由常微分方程的理论可知，只要在区域内连续，且关于满足林普希兹条件，则方程的解存在且唯一。
初值问题的数值解法通常采取“步进法”，而“步进法”又可分为“单步法”和“多步法”两类。
（1）单步法。
所谓“单步法”是指在计算时，只用到前一步的有关信息。
其一般形式为：，主要包括下面三种方法：Euler方法，改进的Euler公式-梯形公式和Runge-Kutta法。
（2）线性多步法。
单步法没有用到前几步计算得到的信息，因此为了提高精度，需重新计算多个点处的函数数值，如RK方法，故计算量较大。
线性多步法的基本思想是充分利用前面的已知信息来构造精度高且计算量小的算法来计算。
多步法常用方法是线性多步法，求解公式为：构造的常用方法是Taylor展开和数值积分方法。
常用的线性多步公式有：四阶Adams显式公式：四阶Adams隐式公式：四阶Milne显式公式：三阶Hamming公式：（隐式公式）预测校正系统和预测校正修正法：一般地，同阶的隐式法比显式法精确，而且数值稳定性好，但隐式公式中的求解较难，需要用到迭代法，这就增加了计算量。
常微分方程初值问题数值解法
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
第五章 常微分方程初值问题数值解法
则有
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章 常微分方程初值问题数值解法
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展 开 法
常微分方程初值问题的数值解法
第七章 常微分方程初值问题的数值解法--------学习小结一、本章学习体会通过本章的学习，我了解了常微分方程初值问题的计算方法，对于解决那些很难求解出解析表达式的，甚至有解析表达式但是解不出具体的值的常微分方程非常有用。
在这一章里求解常微分方程的基本思想是将初值问题进行离散化，然后进行迭代求解。
在这里将初值问题离散化的方法有三种，分别是差商代替导数的方法、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乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
常微分方程初值问题的数值解法
就得到初值问题(7.1),(7.2)的解y(t )的解析表达式。然而
在实际问题和科学研究中所遇到的微分方程往往很复
杂，很多情况下不可能求出它的解析解。有时侯即使
能求出解析解，也会由于很难从解析解中计算函数y(t )
的值而不实用。
例如，容易求出初值问题
y' 1 y cos t,0 t T
注意：这是“折
线法”而非“切
yN
线法”除第一个
点是曲线切线外，
其他点不是切线
而是折线(如右 图所示)。
y2 yy10
t0 a t1 t2
tN b x
§7.2.1 显式单步法的一般形式 显式单步法的一般形式是
yn1 yn h (tn , yn , h), n 0,1, , M 1 (7.2.4)
普希兹)条件。常数L称为函数f 在D0中的Lipschitz常数。
例1 函数f (t, y) t y 在区域D0 (t, y) | 1 t 2, 3 y 4
关于y满足Lipschitz条件，相应的Lipschitz常数可取为L 2
3 存在性定理 定理1 设函数f (t, y)在凸集D R2中有定义，若存在常数
(7.4)
若k 2,则数值解法(7.3)统称为多步法，或具体称 为k步法。
显示法、隐式法与截断误差 若在差分方程(7.3)中，ynk能表示为tn , yn, yn1, , ynk-1, h 的显函数，即
ynk G(tn , yn , yn1, , ynk-1, h), n 0,1, , M - k (7.5)
y0 yn1
Hale Waihona Puke y(t0 yn)
计算方法课件第八章常微分方程初值问题的数值解法
整体截断误差与局部截断误差的关系
定理：如果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 )
常微分方程初值问题数值解法
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) 在一系列离散节点
第7讲 常微分方程初值问题数值解法
yn1 yn hf ( xn1 , yn1 ) n 0,1, 2,, N 1
称上述公式为向后Euler 公式。 向后Euler 公式为隐式格式，需要利用迭代法求解
10
数值分析
第7讲 常微分方程初值问题数值解法
Euler方法的几何意义
③ 图形解
y
•龙格-库塔法
o
x
2
数值分析
第7讲 常微分方程初值问题数值解法
§1 引言 初值问题及其数值解的概念
dy f ( x , y ); a x b 一阶常微分方程初值问题： d x ( ) y( x ) y 0 0
常用的一些 解析解法：
分离变量法、变量代换、 常数变易法、Lapalace变换等
(3)
用Taylor多项式近似并可估计误差
h
2
y( x n 1 ) y( x n h) y( x n ) hy' ( x n )
y' ' ( ) h
2
2! y' ' ( x n )
y( x n ) hy' ( x n )
2!
进一步: 令
y n 1 y( x n 1 ) , y n y ( x n )
计算结果
（1）步长h=0.1的数值解比较表
x 0 0.1
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
精确解 1 1.0048
1.0187 1.0408 1.0703 1.1065 1.1488 1.1966 1.2493 1.3066 1.3679
向前欧拉 1 1
常微分方程初值问题解法
8
4 后退的欧拉方法
(5)
9
(6)
(6)式称为后退的欧拉方法,它是隐式的, 欧拉公式(2)是显式的,
10
(7)
11
12
后退的欧拉方法的局部截断误差:
13
5 梯形方法
(8)
(8)式称为梯形方法.
14
梯形方法的局部截断误差:
15
6. 改进欧拉法及局部截断误差
（1）改进的欧拉公式:
预测步
~ yn 1 yn hf xn , yn
33
提高Runge-Kutta方法的精度的方法
提高积分方法的精度,我们最熟悉的(不一定是最好的)措施是
Richardson 外推法
我们用一个例子予以说明如下 的近似解： Euler法 yn 1 yn hyn y ( h ) ( x) y ( x) c1h c2 h 2 h 将步长减半为 时，有 2 h ( ) 1 1 y 2 ( x) y ( x) c1h c2 h 2 2 4
h) p 1 Yn( (1 1 ch
d
(h) n 1
因此可以从两次计算当中估计出每一步的截断误差,有了这个误差估计之后, 通过与控制误差限比较,就可以控制步长. 注意这个方法增加了计算量.
35
1 1 2 p
) (h (h) 2 Y Y n 1 ; n 1
30
yx 例 求解初值问题ODE : dy dx
易知其精确解为：y 2 2 x x 2 e x
2
,
y (0) 1.
步长都取为 h 0.1 分别用二阶、四阶 R K方法求解：
x
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（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)。
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此前进做出
yp
yn
h( yn
2 xn yn
),
yc
yn
h( y p
2 xn1 ), yp
yn1
1 2
(yp
yc ).
仍取 h ，0.计1 算结果见表9-2. 同例1中欧拉法的计算结 果比较，改进欧拉法明显改善了精度.
24
表9 2 计算结果对比
xn
yn
y(xn ) xn
yn
y(xn )
0.1 1.0959 1.0954 0.6 1.4860 1.4832
（2.12） 则称方法（2.10）具有p阶精度.
若将（2.12）展开式写成
Tn1 (xn , y(xn ))h p1 O(h p2 ), 则 (xn , y(x称n ))为h局p1部截断误差主项.
以上定义对隐式单步法（2.9）也是适用的.
19
对后退欧拉法（2.5）其局部截断误差为
Tn1 y( xn1) y( xn ) hf ( xn1, y( xn1))
)
f ( xn1, yn1 )
hL
y(k) n 1
yn1
.
由此可知，只要 hL 迭1代法（2.6）就收敛到解 .yn1 关于后退欧拉方法的公式误差，从积分公式看到它与欧拉
法是相似的.
12
9.2.2 梯形方法
为得到比欧拉法精度高的计算公式，在等式（2.4）右
端积分中若用梯形求积公式近似，并用 y代n 替 y(xn ), yn1 代替 y(x，n1则) 得
于是可得欧拉法（2.1）的公式误差
y(xn1)
yn1
h2 2
y(n )
h2 2
y(xn ),
称为此方法的局部截断误差.
（2.3）
9
如果对方程（1.1）从 xn到 x积n1分，得
y(xn1) y(xn )
xn1 f (t, y(t))dt.
xn
（2.4）
右端积分用左矩形公式 h f (xn , y近(x似n ，)) 再以 代替 yn y(xn ),代yn替1 也y得(x到n1（) 2.1），局部截断误差也
还可以通过几何直观来考察欧拉方法的精度.
假设 yn ，y即(x顶n )点 落在积Pn分曲线 上，那么，按欧拉方法做出的折线 Pn P便n1是 过点 P的n 切线（图9-2）.
y y(x) y y(x)
图9-2
8
从图形上看，这样定出的顶点 P显n著1 地偏离了原来 的积分曲线，可见欧拉方法是相当粗糙的.
为了控制计算量，通常只迭代一两次就转入下一步的 计算，这就简化了算法.
具体地，先用欧拉公式求得一个初步的近似值 y, n1 称之为预测值，预测值 yn的1精度可能很差，再用梯形公 式（2.7）将它校正一次，即按（2.8）式迭代一次得，这个 结果称校正值，而这样建立的预测-校正系统通常称为改进 的欧拉公式：
Tn1 y(xn1) y(xn ) h(xn , y(xn ), h)
（2.11）
为显式单步法（2.10）的局部截断误差. Tn1 之所以称为局部的，是假设在xn 前各步没有误差.
当 yn y时(x，n )计算一步，则有
y(xn1) yn1 y(xn1) [ yn h (xn , yn , h)] y(xn1) y(xn ) h (xn , y(xn ), h) Tn1.
为了分析计算公式的精度，通常可用泰勒展开将 y(xn1) 在 x处n 展开，则有
y(xn1) y(xn h)
y(xn )
y(xn )h
h2 2
y(n )
n ( xn , xn1).
在 yn y的(x前n )提下，
f (xn , yn ) f (x. n , y(xn )) y(xn )
yn1
yn
h 2
[
f
( xn
,
yn
)
f ( .
梯形方法是隐式单步法，可用迭代法求解.
同后退的欧拉方法一样，仍用欧拉方法提供迭代初值， 则梯形法的迭代公式为
13
y(0) n 1
yn
h
y (k 1) n 1
yn
f (xn , h[f ( 2
yn ); xn , yn
14
如果选取 h充分小，使得
hL 1, 2
则当 k 时有 收敛的.
yn(k1)，这说yn明1 迭代过程（2.8）是
15
9.2.3 单步法的局部截断误差与阶
初值问题(1.1)，(1.2)的单步法可用一般形式表示为
yn1 yn h(xn , yn , yn1, h),
（2.9）
其中多元函数 与 f (x有, y关)，当 含有 时，yn方1法
第9章 常微分方程初值问题数值解法
9.1 引 言
本章要着重考察的一阶方程的初值问题
y f (x, y),
y(
x0
)
y0 .
（1.1） （1.2）
只要函数 f (适x,当y)光滑——譬如关于 满足利y 普 希茨(Lipschitz）条件
f (x, y) f (x, y) L y y .
（1.3）
理论上就可以保证初值问题（1.1），（1.2）的解y y(x) 存在并且唯一.
1
所谓数值解法，就是寻求解 y在(x一) 系列离散节点
x1 x2 xn xn1
上的近似值 y1, y2 , , yn., y相n邻1,两 个节点的间距 hn xn1称为x步n 长.
如不特别说明，总是假定 hi h(i为定1,2数,，) 这时节点为 xn x0 nh(i . 0,1,2, )
h2 2
y(xn ) O(h3 ),
即为（2.3）的结果.
这里 h2
2
y(称xn为) 局部截断误差主
项. 显然 Tn1 O. (h2 )
18
定义2 设 y是(x初) 值问题(1.1)，(1.2)的准确解，若 存在最大整数 p使显式单步法(2.10)的局部截断误差满足
Tn1 y(x h) y(x) h(x, y, h) O(h p1),
yn1
yn
h( yn
2 xn yn
).
取步长 h 0，.1计算结果见表9-1.
表9 1 计算结果对比
xn
yn
y(xn ) xn
yn
y(xn )
0.1 1.1000 1.0954 0.6 1.5090 1.4832
0.2 1.1918 1.1832 0.7 1.5803 1.5492
0.3 1.2774 1.2649 0.8 1.6498 1.6125
22
预测 校正
yn1 yn hf (xn , yn ),
yn1
yn
h 2
[
f
( xn
,
yn
)
f ( xn1, yn1)].
（2.13）
或表为下列平均化形式
y
p
yn
hf
(xn ,
yn ),
yc
yn
hf
( xn1,
y p ),
yn1
1 2
(yp
yc ).
23
例2 用改进的欧拉方法求解初值问题（2.2）. 解 改进的欧拉公式为
初值问题（1.1），（1.2）的数值解法的基本特点是 采取“步进式”，即求解过程顺着节点排列的次序一步一步 地向前推进.
描述这类算法，只要给出用已知信息 yn , yn1, yn2 , 计算 yn的1 递推公式.
首先，要对方程（1.1）离散化，建立求数值解的递
2
推公式. 一类是计算 y时n只1 用到前一点的值 ，称y为n 单步
是（2.3）.
如果在（2.4）中右端积分用右矩形公式 h f (xn1, y(xn1)) 近似，则得另一个公式
yn1 yn hf (xn1, yn1),
（2.5）
称为后退的欧拉法.
后退的欧拉公式与欧拉公式有着本质的区别，后者是 关于 yn的1 一个直接的计算公式，这类公式称作是显式的；