常微分方程初值问题数值解法.
常微分方程初值问题数值解法
常微分方程初值问题的数值解法在自然科学、工程技术、经济和医学等领域中，常常会遇到一阶常微分方程初值问题：(,),,(),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) 该格式是显式，也是二阶方法。
常微分方程初值问题的基本数值解法分析
Ｋ ｅ ｏ ｄｓｏｄｎ ｒ ｉ ｅｅ ｔ ｌ ｑ ａｉｎ ； ｕ ｒｃｌ ｏｕｉｎ；ｏ ｖ ｒｅ ｃ ；ｔｂｌｙ ｙ ｗ ｒ ：ｒ ｉａ ｄｆ ｒｎｉ ｕ ｔ ｓ ｎ ｍｅｉａ ｌｔ ｙ ｆ ａｅ ０ ｓ ｏ ｃ ｎ ｅｇ ｎ ｅ ｓａ ｉｔ ｉ
科学技术及实际生产实践 中的许 多问题都可 归结为微分方程的求解问题 ， 使用较多的是常微分 方程初值 问题的求解． 对于一阶常微分方程的初值 问题 ｄ／ ＝ （，）ｙ ｘ）ｙ，。 ６ 其 中厂 已知 ｙ ｘｆｘ Ｙ ，（。＝ｏ ＜ ， ｄ 为 函数 ，ｏ 初 始值 ． 果 函数 厂 于 变 量 满 足 Ｌｐ ｙ是 如 关 ｉ． ｓｈ ｚ ｃｉ 条件 ， ｔ 则初值 问题有唯一解． 只有当厂 是一些 特 殊 类 型 的 函数 时 ， 能求 出 问题 的解 析 解 ， 一 才 但 般情况下都满足不了生产实践 与科学技术发展的 需要 ， 因此 通常求 其数 值解 法 ．
ｍｅｈ ｄ ａ ｄ ｕ ｄ ｔ ｒ ｎ ｄ ｃ ｅ ｆ ｉｎ ｓ ｍｅｈ ｄ ｔ ｏ ｎ ｎ ｅｅ ｍｉ ｅ ｏ ｆ ｃｅ ｔ ｉ ｔ ｏ ．Ａｔｔ ｅ ｓ ｍｅ ｔ ， ｈ ａ ｅ ｅ ｕ ｅ ｈ ｕ ｅ ｅ ｅ ｏ ｕ ａ ａ ｄ ｔ ｅ ｈ ａ ｉ ｍｅ ｔ ｅ ｐ ｐ ｒｄ ｄ ｃ ｓ ｔ ｅ Ｅ ｌｒ ｓ ｒ ｓ ｆ ｒ ｌ ｎ ｈ ｉ ｍ
ａ ｄ ｍａ ｙ ｏ ｈ ｒｆ ｌ ｓ Ｔｈ ｓｐ ｐ ｒａ ａｙ ｅ ｒ ｅ ｋ ｎ ｓｏ ａ ｉ ｔ ｏ ｓ ｆｒｃ ｎ ｔｕ ｔ ｇ ｎ ｍｅ ｃ ｌｓ ｌ ｔ ｎ ｒｉ ｉ ａ ｎ ｎ ｔ ｅ ｅｄ ． ｉ ａ ｅ ｎ ｚ ｄ ｔ ｅ ｉ ｄ ｆｂ ｓｃ ｍｅｈ ｄ ｏ ｏ ｓｒ ｃ ｉ ｕ ｒ ａ ｏ ｕ ｉ ｓｆ ｎ ｔ ｌ ｉ ｌ ｈ ｎ ｉ ｏ ｏ ｉ
浅谈常微分方程初值问题数值解法
浅谈常微分方程初值问题数值解法在自然科学、工程技术、甚至社会科学的一些领域中，常常会遇见一阶常微分方程的求解问题：()上述问题，寻求解的具体表达式十分困难，仅对一些特殊形式的才有可能找到解的解析表达式，在大多情况下，初值问题的解不能用初等函数表示出来即使可写出解的解析表达式，但因为这些表达式过于复杂，要计算它在某些点上的函数值也异常困难。
在实际问题中，经常需要的恰是解在某些点上的函数值，因此研究初值问题的数值解法十分必要。
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章 常微分方程初值问题数值解法
数值分析
第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 的切线（存在！），则斜率
常微分方程初值问题的数值解法中三种算法的比较
常微分方程初值问题的数值解法中三种算法的比较
常微分方程初值问题的数值解法是数学分析中的一个重要的研究内容，众多的
算法都有助于我们更好地求解一般的初值问题，在这里我们将介绍常微分方程初值问题的三种基本算法，它们是欧拉法、改进欧拉法以及四阶龙格-库塔法。
欧拉法是常微分方程初值问题中最常用的算法，他是一种简洁而又灵活的方法，其基本思想是根据给定的常微分方程和初值，通过积分形式求解精确解，此方法解决的问题比较简单，但它的误差公式与时间步长的N次方有关，误差较大，而且容易出现严重的误差误差，当时间步长To增大时会出现误差振荡。
改进欧拉法是弥补欧拉法缺陷的一种优化算法，它使用线性插值，代替欧拉法
用积分形式计算出来的结果，从而更准确地求出结果，且误差降低，由于它对动态系统的误差有一定的抑制，使得它的运算误差相对于欧拉法是高准确度的，但在某些特殊情况下仍然可能出现误差波动的情况。
四阶龙格-库塔法是在现实生活中最常用的数值解法。
它把问题分解成5种不
同形式的积分公式，并分别交由5个层次的方法来解决，仔细把握每一步的运算，把数值舍入后再运算，虽然该法运算量大，但它的准确性更高，误差相对于其它两种方法要小得多，且具有良好的精度稳定性，具有很好的鲁棒性和适应性，可以很好地用于对解初值问题作出估计和预测。
综上，这三种数值解法都有自身的特点，欧拉法计算简单，但误差较大；改进
欧拉法的精度和误差抑制能力更强；四阶龙格-库塔法的算术精度更高，出现误差
波动的概率最低，在可靠性方面更加准确。
因此，应用的时机对于三种算法的选择就显得尤为重要。
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解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乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
数值分析李庆扬第9章常微分方程初值问题数值解法讲义.
② 由 x0 , y0 f x0 , y0 切线 P0P1 ，
切线与 x x1 交点 P1 ： y1 的近似值 ；
③ 再由 x1 , y1 向前推进到 P2 ， 得到折线 P0P1 Pn ，近似 y yx 。
7
2021年5月4日
《数值分析》 黄龙主讲
h
yxn
yxn1
yn1 yn f
h
xn1 , yn1
yn1 yn h f xn1 , yn1
——后退的欧拉公式（隐式）
注意：① 显式计算方便，隐式稳定性较好；
② 上式隐含 yn1 ? ，采用迭代法求解。
12
2021年5月4日
《数值分析》 黄龙主讲来自欧拉公式的另一种理解：
将常微分方程 y f x, y 改写 dy f t , ytdt
“步进式”：顺着节点排列顺序，一步一步地向前推进。
步长：常用等步长 hn xn1 xn ，节点为 xn x0 nh 单步法：计算 yn1 时，只用到前一点的值 yn k 步法：计算 yn1 时，用到前面 k 点的值 yn , yn1 , , ynk1
5
2021年5月4日
《数值分析》 黄龙主讲
对微分方程从 xn 到 xn1 积分
y xn1 yxn
xn1 f t , yt dt
xn
由积分左矩形公式得
xn1 xn
f
t ,
yt dt
hf
xn ,
yxn
例如：
lim
h0
yxn1
h
yxn
yxn
yxn1
h
yxn
yxn
f xn , yxn
常微分方程初值问题的数值积分法
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.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.改进
常微分方程初值问题的数值解法
--------学习小结一、本章学习体会通过本章的学习，我了解了常微分方程初值问题的计算方法，对于解决那些很难求解出解析表达式的，甚至有解析表达式但是解不出具体的值的常微分方程非常有用。
在这一章里求解常微分方程的基本思想是将初值问题进行离散化，然后进行迭代求解。
在这里将初值问题离散化的方法有三种，分别是差商代替导数的方法、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 T y 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 ε+=。
常微分方程初值问题的的数值解法
本章讨论常微分方程初值问题的数值解法
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) 在一系列离散节点
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（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}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。
实验八 常微分方程初值问题数值解法报告
实验八 常微分方程初值问题数值解法一、基本题科学计算中经常遇到微分方程（组）初值问题，需要利用Euler 法，改进Euler 法，Rung-Kutta 方法求其数值解，诸如以下问题：（1） ()⎪⎩⎪⎨⎧=-='004y xy y x y 20≤<x分别取h=0.1,0.2,0.4时数值解。
初值问题的精确解245x y e -=+。
（2） ()⎩⎨⎧=--='0122y y x y 01≤≤-x用r=3的Adams 显式和预 - 校式求解取步长h=0.1，用四阶标准R-K 方法求值。
（3）()()()100010321331221==-='⎪⎩⎪⎨⎧-='-='='y y y y y y y y y 10≤≤x用改进Euler 法或四阶标准R-K 方法求解取步长0.01,计算(0.05),(0.1y y y 数值解，参考结果 123(0.15)0.9880787,(0.15)0.1493359,(0.15)0.8613125y y y ≈-≈≈。
（4）利用四阶标准R- K 方法求二阶方程初值问题的数值解（I ）()()⎩⎨⎧='==+'-''10,00023y y y y y 02.0,10=≤≤h x(II)()()()⎩⎨⎧='==+'--''00,10011.02y y y y y y 1.0,10=≤≤h x(III)()()⎪⎩⎪⎨⎧='=+='00,101y y e y y x 1.0,20=≤≤h x(IV)()()⎩⎨⎧='==+''00,100sin y y y y 2.0,40=≤≤h x二、应用题1. 小型火箭初始质量为900千克，其中包括600千克燃料。
火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
常微分方程初值问题数值解法朱欲辉（浙江海洋学院数理信息学院, 浙江舟山316004）[摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点.[关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计Numerical Method for Initial-Value ProblemsZhu Yuhui(School of Mathematics, Physics, and Information Science,Zhejiang Ocean University, Zhoushan, Zhejiang 316004)[Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis.[Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate1前言自然界和工程技术中的很多现象, 例如自动控制系统的运行、电力系统的运行、飞行器的运动、化学反应的过程、生态平衡的某些问题等, 都可以抽象成为一个常微分方程初值问题. 其真解通常难以通过解析的方法来获得, 至今有许多类型的微分方程还不能给出解的解析表达式, 一般只能用数值的方法进行计算. 有关这一问题的研究早在十八世纪就已经开始了, 特别是计算机的普遍应用, 许多微分方程问题都获得了数值解, 从而能使人们认识解的种种性质及其数值特征, 为工程技术等实际问题提供了定量的依据.关于常微分方程初值问题的数值计算方法, 许多学者己经做了大量的工作. 1768年, Euler 提出了关于常微分方程初值问题的方法, 1840年, Cauchy 第一次对初值问题进行了仔细的分析, 早期的常微分方程数值解的问题来源于天体力学. 在1846年, 当Adams 还是一个学生的时候, 和Le Verrier 一起根据天王星轨道中出现的己知位置, 预测了它下一次出现的位置. 1883年, Adams 提出了Adams 一Bashforth 和Adams 一Moulton 方法. Rull (1895年)、Heun(1900年)和Kutta (1901年)提出Runge.Kutta 方法.二十世纪五十年代, Dahlquist 建立了常微分方程数值解法的稳定性理论, 线性多步法是常微分方程初值问题的一种数值方法. 由于通常的数值方法, 其绝对稳定区域是有限的, 不适用于求解刚性常微分的初值问题. 刚性微分方程常常出现于航空、航天、热核反应、自动控制、电子网络及化学动力学等一系列与国防和现代化建设密切相关的高科技领域, 具有无容置疑的重要性. 因此, 刚性微分方程的研究工作早在二十世纪五十年代就开始了, 1965年, 在爱丁堡举行的IFIP 会议后, 更进一步地认识刚性方程的普遍性和重要性. 自从六十年代初, 许多数值分析家致力于探讨刚性问题的数值方法及其理论, 注意到刚性问题对传统数值积分方法所带来的挑战. 这一时期, 人们的研究主要集中在算法的线性稳定性上, 就是基于试验方程y y C λλ=∈，，（）数值解的稳定性研究. 在此领域发表了大量的论文, 取得了许多重要的理论成果. 例如, 1963年, Dahlquist 给出A 稳定性理论, 1967年, Widlund 给出()A α—稳定性理论, 1969年, Gear 将A —稳定性减弱, 给出刚性(Stiff )稳定性理论, 并找到了当k 6≤的k 步k 阶的刚性稳定方法, 1969年Dill 找到刚性稳定的7阶和8阶以及1970年Jain 找到刚性稳定的9阶到11阶, 但可用性没有检验. 这些稳定性理论和概念都是在线性试验方程的框架下推导出的, 从严格的数学意义上来说, 这些理论只适用于常系数线性自治系统. 但从实用的观点来说, 这些理论无疑是合理和必要的, 对刚性问题的算法设计具有重要的指导意义. 在八十至九十年代, 国内也有一些学者研究线性理论, 主要有匡蛟勋、陈果良、项家祥、李寿佛、黄乘明、李庆扬和费景高等.线性理论虽然对一般问题具有指导作用, 但其不能作为非线性刚性问题算法的稳定性理论研究基础. 为了将线性理论推广到非线性问题中, 人们开始对非线性模型问题进行研究.但是, 早期文献主要致力于数值方法基于经典Lipschitz条件下的经典收敛理论, 即认为良好的稳定性加上经典相容性和经典相容阶就足以描述方法的整体误差性态. 直到1974年, Prothero和Robinson首先注意到算法的经典误差估计由于受刚性问题巨大参数的影响而严重失真, 产生阶降低现象, 这时人们认识到经典收敛理论对于非线性刚性问题以及线性模型的不足. 于是, 1975年, Dahlquist和Butcher分别提出了单支方法和线性多步法的G一稳定概念和B稳定概念. 这两个概念填补了非线性稳定性分析理论, 引起了计算数学家们的极大关注, 在上述理论的基础上, 1975年至1979年, Burrage和Butcher提出了AN一稳定性与BN一稳定性概念, 并相应地建立了基本的B一稳定及代数稳定理论. 1981至1985年, Frank, schneid和ueberhuber建立Runge一kutta方法的B一收敛理论. B稳定与B一收敛理论统称B一理论, 它是常微分方程数值解法研究领域的巨大成就之一, 是刚性问题算法理论的突破性进展, 标志着刚性问题研究从线性向非线性情形深入发展. 国内也有众多学者致力于B一理论的研究, 如李寿佛、曹学年等. 1989年, 李寿佛将Dahlquist的G一稳定概念推广到更一般的(C,P,Q)代数稳定, 克服了G稳定的线性多步法不能超过二阶的限制. 对于一般线性方法, 李寿佛建立了一般线性方法的(K,P,Q)稳定性理论及(K,P,Q)弱代数稳定准则和多步Runge-Kutta法的一系列代数准则. 此外, Dahquist, Butcher和Hairer分别深刻地揭示了单支方法、一般线性方法和Runge.Kutta方法线性与非线性稳定性之间的内在关系. 为了求解刚性微分方程, 不少文献中构造含有稳定参数的线性多步方法, 利用适当选择稳定参数来扩大方法的稳定区域. 所有改进的思想, 都是通过构造一些特殊的显式或隐式线性多步法,Aα一稳定的α角增大. 八十年代, 就成为国内外学者所使其具有增大的稳定域, 或使()研究的一个课题, 学者主要有Rodabaugh和Thompson、Feinberg、李旺尧、李寿佛、包雪松、徐洪义、刘发旺、匡蛟勋、项家祥、蒋立新、李庆扬、谢敬东和李林忠等.当前国内外研究刚性问题的一个主要趋势就是在B一理论指导下寻找更为有效的新算法. 另一个发展趋势就是力图突破单边lipschitz常数和内积范数的局限, 建立比B一理论更为普遍的定量分析收敛理论. 近年来, 刚性延迟系统的算法研究成为刚性问题的另一个热点研究领域, 张诚坚将Burrage等人创立的针对刚性常微分系统的B一理论拓展到非线性刚性延迟系统.常微分方程的数值算法发展到今天己有了线性多步法、龙格一库塔法和在此基础上发展起来的单支方法、分块方法、循环方法、外推法、混合方法、二阶导数法以及各种常用的估校正算法. 其中经常用到的线性多步法公式有Euler 公式、Heun 公式、中点公式、Milne 公式、Adams 公式、simpson 公式、Hamming 公式, Gear 方法、Adams 预估一校正法和Mile 预估一Hamming 校正法公式等, 此外还包含许多迄今尚末探明的新公式. Burage 曾将线性多步法和Runge —Kutta 法比作大海中的两座小岛, 在浩瀚的汪洋之中, 还有许多到现在没有发现的新方法.本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler 法、改进的欧拉法、显式龙格-库塔法、隐式龙格-库塔法以及线性多步法中的Adams 显隐式公式和预测校正公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 进行了计算机程序算法的分析与实现, 以计算机的速度优势来弥补计算量大的不足. 从得出的结果对比各种方法的优缺点.2常微分方程初值问题的数值解法2.1 数值方法的基本思想考虑一阶常微分方程初值问题(,),[,],dx f x y x a b dy=∈ 00()y x y = (2.1) 的数值解法, 其中f 是x 和y 的己知函数, 0y 是给定的初始值.对于常微分方程初值问题(2.1)数值解法, 就是要算出精确解()y x 在区间[,]a b 上的一系列离散节点011n n a x x x x b -=<<<<=的函数值0()y x , 1()y x , ,()n y x 的近似值0y , 1y , , n y .相邻两个节点的间距1i i h x x +=-称为步长, 这时节点也可以表示为0i x x ih =+(1,2,,)i n =. 数值解法需要把连续性的问题加以离散化, 从而求出离散节点的数值解. 通常微分方程初值问题(2.1)的数值方法可以分为两类:(1) 单步法-----计算()y x 在1n x x +=处的值仅取决于n x 处的应变量及其导数值.(2) 多步法-----计算()y x 在1n x x +=处的值需要应变量及其导数在1n x +之前的多个网个节点出的值.2.2 Euler 方法2.2.1 Euler 公式若将函数()y x 在点n x 处的导数'()n y x 用两点式代替, 即'1()()()n n n y x y x y x h+-≈, 再用n y 近似地代替()n y x , 则初值问题(2.1)变为100(,)(),0,1,2,n n n n y y hf x y y y x n +=+⎧⎨==⎩ (2.2) (2.2)式就是著名的欧拉(Euler)公式.2.2.2 梯形公式欧拉法形式简单精度低, 为了提高精度, 对方程'(,)y f x y =的两端在区间1[,]n n x x +上积分得11()()[,()]n n x n n x y x y x f x y x dx ++=+⎰(2.3)改用梯形方法计算其积分项, 即 1111[,()][(,())(,())]2n n x n n n n n n x x x f x y x dx f x y x f x y x ++++-≈+⎰ 代入式(2.3), 并用n y 近似代替式中()n y x 即可得到梯形公式111[(,)(,)]2n n n n n n h y y f x y f x y +++=++ (2.4) 由于数值积分的梯形公式比矩形公式精度高, 所以梯形公式(2.4)比欧拉公式(2.2)的精度高一个数值方法.式(2.4)的右端含有未知的1n y +, 它是关于1n y +的函数方程, 这类方法称为隐式方法.2.2.3 改进的欧拉公式梯形公式实际计算时要进行多次迭代, 因而计算量较大. 在实用上, 对于梯形公式(2.4)只迭代一次, 即先用欧拉公式算出1n y +的预估值1n y +, 再用梯形公式(2.4)进行一次迭代得到校正值1n y +, 即1111(,)[(,)(,)]2n n n n n n n n n n y y hf x y h y y f x y f x y ++++⎧=+⎪⎨=++⎪⎩ (2.5) 2.2.4 欧拉法的局部截断误差衡量求解公式好坏的一个主要标准是求解公式的精度, 因此引入局部截断误差和阶数概念.定义 2.1 在n y 准确的前提下, 即()n n y y x =时, 用数值方法计算1n y +的误差111()n n n R y x y +++=-, 称为该数值方法计算1n y +时的局部截断误差.对于欧拉公式, 假定()n n y y x =, 则有'1()[(,())]()()n n n n n n y y x h f x y x y x hy x +=+=+而将真解()y x 在n x 处按二阶泰勒展开式有2''''11()()(),(,)2!n n n n n h y y x hy x y x x ξξ++=++∈ 因此有2''11()()2!n n h y x y y ξ++-= 定义 2.2 若数值方法的局部截断误差为1()p O h +, 则称这种数值方法的阶数是P. 步长(h<1)越小, P 越高, 则局部截断误差越小, 计算精度越高.2.3 Runge －Kutta 方法2.3.1 Runge －Kutta 方法的基本思想欧拉公式可改写成111(,)n n n n y y hK K f x y +=+⎧⎨=⎩则1n y +的表达式与1()n y x +的泰勒展开式的前两项完全相同, 即局部截断误差为2()O h .改进的欧拉公式又可改写成 1121211()2(,)(,)n n n n n n h y y K K K f x y K f x y hk ++⎧=++⎪⎪=⎨⎪=+⎪⎩. 上述两组公式在形式上有一个共同点: 都是用(,)f x y 在某些点上值的线形组合得出1()n y x +的近似值1n y +, 而且增加计算(,)f x y 的次数, 可提高截断误差的阶. 如欧拉公式, 每步计算一次(,)f x y 的值, 为一阶方法. 改进的欧拉公式需计算两次(,)f x y 的值, 它是二阶方法. 它的局部截断误差为3()O h .于是可考虑用函数(,)f x y 在若干点上的函数值的线形组合来构造近似公式, 构造时要求近似公式在(,)n n x y 处的泰勒展开式与解()y x 在n x 处的泰勒展开式的前面几项重合, 从而使近似公式达到所需要的阶数. 既避免求偏导, 又提高了计算方法精度的阶数. 或者说, 在1[,]n n x x +这一步内多预报几个点的斜率值, 然后将其加权平均作为平均斜率. 则可构造出更高精度的计算格式, 这就是Runge －Kutta 方法的基本思想. 2.3.2 Runge －Kutta 方法的构造一般地, Runge －Kutta 方法设近似公式为 （下面的公式修改了）11111(,)(,)p n n i i i n n i i n i n ij j j y y h c K K f x y K f x a h y h b K +=-=⎧=+⎪⎪⎪=⎨⎪⎪=++⎪⎩∑∑ ()2,3,,i p =其中i a , ij b , i c 都是参数, 确定它们的原则是使近似公式在m y 处的泰勒展开式与()y x 在n x 处的泰勒展开式的前面的项尽可能多的重合, 这样就使近似公式有尽可能高的精度. 以此我们可以通过一个复杂的计算过程得出常用的的三阶和四阶Runge －Kutta 公式1123121312(4)6(,)(,)22(,2)n n n n n n n n h y y K K K K f x y h h K f x y K K f x h y hK hK +⎧=+++⎪⎪=⎪⎨⎪=++⎪⎪=+-+⎩ 和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 +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩(2.6) 式(2.6)称为经典Runge －Kutta 方法.2.4 线性多步法在逐步推进的求解过程中, 计算1n y +之前事实上已经求出了一系列的近似值0y , 1y, , n y . 如果充分利用前面多步的信息来预测1n y +, 则可期望获得较高的精度, 这就是构造多步法的基本思想.线性k 步方法的一般公式为11101(,)k k n j n j j n j n j j j y a y h b f x y --+---==-=+∑∑ (2.7)其中j a , j b 均为与n 无关的常数, 110k k a b --+≠. 当10b -=时为显格式; 当10b -≠时为隐格式. 特别当0011,1,0k a b b -====时为Euler 公式；当00111,1,2k a b b -====时为梯形公式. 定义 2.3 称1111101()[(,)]k k n n n j n j j n j n j j j R y x y a y h b f x y --+++---==-=-=+∑∑为k 步公式(2.7)在1n x +处的局部截断误差. 当11()p n R O h ++=时称式(2.7)是p 阶的.应用方程'()(,())y x f x y x =可知局部截断误差也可写成为 11'11101()[()]k k n n n j n j j n j j j R y x y a y h b y x --+++--==-=-=+∑∑定义2.4 如果线性k 步方法(2.7)至少是1阶的, 则称是相容的; 如果线性k 步法(2.7)是p 阶的, 则称是p 阶相容的.2.4.1 Adams 外插法将微分方程(,),[,],dx f x y x a b dy=∈ 的两端从n x 到1n x +进行积分, 得到11()()(,())n n x n i x y x y x f x y x dx ++=+⎰(2.8)我们用插值多项式代替右端的被积函数. Adams 外插法选取k 个点1,1,,,n n n k x x x -=-+作为插值基点构造(,)f x y 的k-1阶多项式Adams 外插法的计算公式为110k j n n j m j y y h a f -+==+∇∑其中j a 满足如下代数递推式:j 120111a 1231j j a a a j --++++=+, 0,1,j =根据此递推公式, 可逐个的计算j a ()0,1,j =, 表2.1给出了j a 的部分数值:表2.12.4.2 Adams 内插法根据插值理论知道, 插值节点的选择直接影响着插值公式的精度, 同样次数的内插公式的精度要比外插公式的高.仍假定已按某种方法求得问题(2.1)的解()y x 在(0,1,)i x i n =处的数值i y , 并选取插值节点1,,,,n k n n p x x x -++, p 是正整数, 用Lagrang e 型插值多项式(),1()p n k L x -构造'(,)y f x y =可以导出解初值问题(2.1)的Adams 内插公式:1()1,1()()()n n x p n n n k x y x y x L x dx ++-=+⎰ (2.9) 当0p =时上式就退化为内插公式.内插公式(2.9)除了包含f 在1,,n k n x x -+处的已知值外, 还包含了在点,,n n p x x +, 处的未知值.因此内插公式(2.9)只给出了未知量,,n n p y y +的关系式, 实际计算时仍需要解联立方程组. 1p =的内插公式是最适用的, (1),1()n k L x -采用Newton 向后插值公式得到Adams 内插公式*10kj n n j m j y y h a f +==+∇∑其中系数*j a 定义为0*1(1)j j a d j ττ--⎛⎫=- ⎪⎝⎭⎰0,1,j=其中*j a 满足如下代数递推式:****1201,01110,0231j j j j a a a a j j --=⎧++++=⎨>+⎩根据此递推公式, 可逐个的计算*ja()0,1,j =, 表2.2给出了*j a 的部分数值:表2.22.5 算法的收敛性和稳定性 2.5.1 相容性初值问题(2.1):(,),[,],dxf x y x a b dy=∈ 00()y x y =的显式单步法的一般形式为1(,,)n n n n y y h x y h ϕ+=+, 0,1,,1n N =- (2.10)其中b ah N-=, n x a nh =+. 这样我们用差分方程初值问题(2.10)的解n y 作为问题(2.1)的解()n y x 在n x x =处的近似值, 即()n n y x y ≈.因此, 只有在问题(2.1)的解使得()()(,(),)y x h y x x y x h hϕ+--逼近'(,())0y f x y x -=时, 才有可能使(2.10)的解逼近于问题（2.1）的解. 从而, 我们期望对任一固定的[,],x a b ∈ 都有()()lim[(,(),)]0h y x h y x x y x h hϕ→+--=假设增量函数(,,)x y h ϕ关于h 连续, 则有(,,)(,)x y h f x y ϕ=定义 2.5 若关系式(,,)(,)x y h f x y ϕ=成立, 则称单步法(2.10)与微分方程初值问题(2.1)相容.容易验证, Euler 法与改进Euler 法均满足相容性条件. 事实上, 对Euler 法, 增量函数为(,,)(,)x y h f x y ϕ=自然满足相容性条件, 改进Euler 法的增量函数为1(,,)[(,)(,(,))]2x y h f x y f x h y hf x y ϕ=+++因为(,)f x y , 从而有1(,,)[(,)(,)](,)2x y h f x y f x y f x y ϕ=+=所以Euler 法与改进Euler 法均与初值问题（2.1）相容. 一般的, 如果显示单步法有p 阶精度(0)p >, 则其局部误差为1()[()(,(),)]()p y x h y x h x y x h O h ϕ++-+=上式两端同除以h , 得10n ε+→令0h →, 如果(,,)x y h ϕ, 则有'()(,,)0y x x y h ϕ-=即(,,)(,)x y h f x y ϕ=所以0p >的显示单步法均与初值问题(2.1)相容.由此各阶RK 方法与初值问题(2.1)相容.2.5.2 收敛性定义2.6 一种数值方法称为是收敛的, 如果对于任意初值0y 及任意(,],x a b ∈ 都有lim ()n h y y x →= ()x a nh =+其中()y x 为初值问题(2.1)的准确解.按定义, 数值方法的收敛性需根据该方法的整体截断误差1n ε+来判定. 已知Euler 方法的整体截断误差有估计式()11()2L b a n hM e O h Lε-+≤-= 当0h →, 10n ε+→, 故Euler 方法收敛.定理 2.1 设显式单步法具有p 阶精度, 其增量函数(,,)x y h ϕ关于y 满足Lipschitz 条件, 问题(2.1)是精确的, 既00()y x y =, 则显式单步法的整体截断误差为111()()p n n n y x y O h ε+++=-=.定理 2.2设增量函数(,,)x y h ϕ在区域S 上连续, 且关于y 满足Lipschitz 条件时, 则显式单步法收敛的充分必要条件是相容性条件成立, 即(,,)(,)x y h f x y ϕ=.2.5.3 稳定性定义 2.7 如果存在正常数0h 及C , 使得对任意的初始出发值00,y y , 单步法(2.10)的相应精确解,n n y y , 对所有的00h h <≤, 恒有00n n y y C y y -≤-, nh b a ≤-则称单步法是稳定的.定理2.3若(,,)x y h ϕ对于[,]x a b ∈, 0(0,]h h ∈以及一切实数y, 关于y 满足Lipschitz 条件, 则单步法(2.7)是稳定的.定义 2.8对给定的微分方程和给定的步长h, 如果由单步法计算n y 时有大小为δ的误差, 即计算得出n n y y δ=+, 而引起其后之值()m y m n >的变化小于δ, 则说该单步法是绝对稳定的.一般只限于典型微分方程'y y μ=考虑数值方法的绝对稳定性, 其中μ为复常数(我们仅限于μ为实数的情形). 若对于所有的(,)h μαβ∈, 单步法都绝对稳定, 则称(,)αβ为绝对稳定区间.根据以上定义我们可以得出Euler 方法的绝对稳定区间为(-2,0), 梯形公式的绝对稳定区间为(,0)-∞, Runge －Kutta 的绝对稳定区间为(-2.78,0).3 实例分析例3.1分别用Euler 法、改进Euler 法、Runge －Kutta 解初值问题'22(0 1.2)(0)1y xy x y ⎧=-≤≤⎨=⎩的数值解, 取h=0.1, N=10, 并与真实解作出比较.由于该简单方程可以用数学方法求得其精确描述式211y x =+, 所以可以据此检验近似数值解同真实解的误差情况. 对于其他一些结构复杂的常微分方程的数值解实现方法也是一样的.（1）使用欧拉算法进行一般求解经过欧拉算法程序计算得出i y 在各点i x 处的近似数值解及各自同精确解的误差, 如下表3.1表3.1从实验结果看误差已经不算太小, 况且这还仅仅是一个用于实验的比较简单的常微分方程, 对于实际工程中应用的结构复杂的方程其求解结果的误差要远比此大得多, 由于还存在着局部截断误差和整体截断误差, 因此可采取一定措施来抑制减少误差, 尽量使结果精确.（2）使用改进欧拉算法进行一般求解经过改进欧拉算法程序计算得出i y 在各点i x 处的近似数值解及各自同精确解的误差,如下表3.2可以看出, 这种经过改进的欧拉算法所存在的误差已大为减少, 误差的减少主要是由于先利用了欧拉公式对1i y +的值进行了预估, 然后又利用梯形公式对预估值作了校正, 从而在预估—校正的过程中减少了误差. 此算法有一定的优点, 在一些实际而且简单的工程计算中可以直接单独应用.（3）使用经典Runge －Kutta 法进行一般求解经过经典Runge －Kutta 法程序计算得出i y 在各点i x 处的近似数值解及各自同精确解的误差, 如下表3.3表3.3从表3.3记录的程序运行结果来看, 所计算出来的常微分方程的数值解是非常精确的, 用浮点数进行运算时由近似所引入的误差几乎不会对计算结果产生影响, 这样高的精度是很令人满意的. 无论是从计算速度上还是从计算精度要求上, 都能很好地满足工程计算的需要.例3.2分别用Euler 法、改进Euler 法、Runge －Kutta 解初值问题'22(0 1.2)(0)1y xy x y ⎧=-≤≤⎨=⎩的数值解, 取等分数分别为N=12, 24, 36, … ,120, 分别计算出与真实解的最大误差.（1） 使用欧拉算法进行一般求解经过计算我们可以得出初值问题在取不通步长时, 数值解与真实解的最大误差如表3.4表3.4从表3.4记录的程序运行结果来看当等分数N 变大时, 它的误差正在减小, 根据定义2.6我们可以证得该方法是收敛的. （2）使用改进欧拉算法进行一般求解经过计算我们可以得出初值问题在取不通步长时, 数值解与真实解的最大误差如表3.5表3.5从表3.5记录的程序运行结果来看当等分数N 变大时, 它的误差正在减小, 根据定义2.6我们可以证得该方法是收敛的.（3）使用经典Runge －Kutta 法进行一般求解经过计算我们可以得出初值问题在取不通步长时, 数值解与真实解的最大误差如表3.6表3.6从表3.6记录的程序运行结果来看当等分数N变大时,它的误差正在减小, 根据定义2.6我们可以证得该方法是收敛的.（4）收敛速度对比10-的最小等分为对比以上三种方法的收敛速度, 我们计算出了它们的误差精度达到5数N如下表3.7表3.7从表3.7记录的程序运行结果来看Runge－Kutta法的收敛速度最快, 改进Euler法其次, 而Euler法最差. 由此看来Runge－Kutta法是他们当中最理想的数值解法.小结针对工程计算中遇到的常微分方程初值问题的求解, 根据实际情况引如了计算机作为辅助计算工具, 并根据高等数学的有关知识将欧拉公式、改进的欧拉公式、经典龙格－库塔方法等融入到计算机程序算法之中, 充分利用了计算机的速度优势, 大大减轻了工程技术人员的劳动强度, 同时也使计算结果更加可靠、准确。