数值分析(研究生)第七章常微分方程的数值解法二.
合集下载
数值分析常微分方程数值解
许多实际问题的数学模型是微分方程或微分方程的定解问题。
如物体运动、电路振荡、化学反映及生物群体的变化等。
常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。
若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。
因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dxdy， (9-1) 的数值解法具有典型性。
常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。
用解析方法只能求出线性常系数等特殊类型的方程的解。
对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。
因此研究微分方程的数值解法是非常必要的。
只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。
由常微分方程的理论知，如果(9-1)中的),(y x f 满足条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在上关于满足Lipschitz 条件，即存在常数，使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。
在下面的讨论中，我们总假定方程满足以上两个条件。
所谓数值解法，就是求问题(9-1)的解)(x y y =在若干点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的方法。
),,2,1(N n y n =称为问题(9-1)的数值解，n n x x h -=+1称为由到1+n x 的步长。
今后如无特别说明，我们总假定步长为常量。
建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数在问题(9-1)中，若用向前差商hx y x y n n )()(1-+代替)(n x y '，则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样，问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到，按式(9-2)由初值经过步迭代，可逐次算出N y y y ,,21。
数值分析常微分方程数值解法
7
第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
第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
数值分析 (77) 第7章 常微分方程数值解法PPT课件
yi+1)中，因此?，通常用迭代法求解。具体做法是：先
《 数
用欧拉公式(6―3)?求出一个y(0)i+1作为初始近似，然后
值 分
再用改进的欧拉公式(6―5)进行迭代，即
析
》
yi(01) yi hf (xi, yi )
yi(k11)
yi
h[ 2
f
(xi,
yi )
f
(xi1,
y(k) i1
)
k 0,1,2,
y0，y1，y2，…，yn
常取离散点x0，x1，x2，…，xn为等距，即
《
x i+1-xi=h，i=0，1，2，…，n-1
数 值
h称为步长。图61表示为初值问题(6―1)在n+1个
分 析
离散点上的准确解y(x)的近似值。
》
第6章 常微分方程数值解法
§2 欧拉法和改进的欧拉法
2.1 欧拉法(折线法)
数 值
法进行改进，将在一点(xi，yi)的切线斜率f(xi?，yi)用
分 两点的平均斜率来代替，即
析
》
f(xi,yi)1 2[f(xi,yi)f(xi 1,yi 1)]
第6章 常微分方程数值解法
代入(6―3)式得
yi1
yi
h[f 2
(xi,
yi)
f
(xi1,
yi1)]
i 0,1,2, ,n1
(6―5)
《 数 值 分 析 》
第6章 常微分方程数值解法
《 数 值 分 析 》
图 6.2
第6章 常微分方程数值解法
《 数 值 分 析 》
图 6.3
第6章 常微分方程数值解法
2.2 改进的欧拉法
第七章常微分方程数值解法
h2 h3 y ( xi 1 ) y ( xi h) y ( xi ) hy '( xi ) y ''( xi ) y '''( xi ) 2! 3!
丢掉高阶项，有
y( xi 1 ) y( xi h) y( xi ) hy '( xi ) yi hf ( xi , yi )
| f ( x, y1 ) f ( x, y2 ) | L | y1 y2 | ,
那么模型问题在 [ a, b] 存在唯一解。
Lipschitz 连续： | f ( x, y1 ) f ( x, y2 ) | L | y1 y2 | .
(1) 比连续性强: y1 y2 可推出 f ( x, y1 ) f ( x, y2 ) ; (2) 比连续的 1 阶导弱：具有连续的 1 阶导，则
f | f ( x, y1 ) f ( x, y2 ) || ( ) || y1 y2 | L | y1 y2 | . y
常微分方程数值解法
目标：计算出解析解 y ( x) 在一系列节点 a x0 x1 xn1 xn b 处的近似值 yi y( xi ) ，即所谓的数值解。节点间距 hi xi 1 xi ，一般 取为等距节点。
常微分方程初值问题的数值解法一般分为两大类： (1)单步法：在计算 yn 1 时，只用到前一步的值，即用到 xn1 , xn , yn ，则给定初
值之后，就可逐步计算。例如 Euler 法、向后欧拉法、梯形公式、龙格-库塔法；
(2) 多步法： 这 类 方 法 在 计算 yn 1 时 ， 除 了 用 到 xn1 , xn , yn 外 ， 还 要 用到
数值分析(研究生)第七章常微分方程的数值解法二
化作一阶微分方程组求解. 引入新变量 y1 = y, y2 = y, ... , yn = y( n-1)
y1 y -1 n n y = y2 . . . = yn = f ( x , y1 , ... , yn )
初值条件为：
y1 ( x0 ) = a0 y2 ( x0 ) = a1 ... y n ( x 0 ) = a n -1
上页 下页
返回
亚当姆斯显式公式
利用k+1 个节点上的被积函数值 f i , f i -1 , ... , f i -k 构造 k 阶牛顿 后插多项式 N k ( xi + t h) , t [0, 1], 有
xi +1 xi
f ( x, y( x))dx = N k ( xi + t h) h dt + Rk ( xi + t h) h dt
§5 收敛性和稳定性
一、收敛性
定义 若某算法对于任意固定的 x = xi = x0 + i h，当 h0
( 同时 i ) 时有 yi y( xi )，则称该算法是收敛的. 例3
y = y 就初值问题 y (0) = y0
考察欧拉显式格式的收敛性.
解：该问题的精确解为 y( x) = y0e x 尤拉公式为 yi +1 = yi + h yi = (1 + h) yi limy1 + 1 + /hh=i e 0 ( i = ( h)1 ) y
y( x ) = f ( x , y ) y ( x 0 ) = y0
前述所有公式皆 适用于向量形式.
上页 下页
返回
常微分方程数值解法5262115页PPT文档
x 1 ( t ) 表示时刻 t 食饵的密度，x 2 ( t ) 表示捕食者的密度；
r 表示食饵独立生存时的增长率；
d 表示捕食者独立生存时的死亡率；
a 表示捕食者的存在对食饵增长的影响系数，反映捕
食者对食饵的捕获能力；
b 表示食饵的存在对捕食者增长的促进系数，反映食
饵对捕食者的喂养能力
150 100
令 y 1 y ,y 2 y ',y 3 y '', ,y n y ( n 1 )
可以将以上高阶微分方程化为如下一阶常微分方程组
y1 ' y2 y2 ' y3 yn ' an(x)y1
a1(x)yn f (x)
例：P120,1(a),Bessel方程
常微分方程的数值解
一般地，凡表示未知函数，未知函数的导 数与自变量之间的关系的方程叫做微分方 程．未知函数是一元函数的，叫常微分方 程;未知函数是多元函数的，叫做偏微分方 程．
如
y ' x y'x2y2 y''y'xy
Matlab实现 [t,x]=ode45(f,ts,x0,options,p1,p2,......)
50 0 0
30 20 10
0 0
10
20
50
30
20
10
0
30
0
10
8
6
4
2
100
0
50
100
150
50
100
高阶常微分方程的解法
高阶常微分方程
y ( n ) a 1 ( x ) y ( n 1 ) a ( n 1 ) ( x ) y ' a n ( x ) y f( x )
r 表示食饵独立生存时的增长率；
d 表示捕食者独立生存时的死亡率；
a 表示捕食者的存在对食饵增长的影响系数，反映捕
食者对食饵的捕获能力；
b 表示食饵的存在对捕食者增长的促进系数，反映食
饵对捕食者的喂养能力
150 100
令 y 1 y ,y 2 y ',y 3 y '', ,y n y ( n 1 )
可以将以上高阶微分方程化为如下一阶常微分方程组
y1 ' y2 y2 ' y3 yn ' an(x)y1
a1(x)yn f (x)
例：P120,1(a),Bessel方程
常微分方程的数值解
一般地，凡表示未知函数，未知函数的导 数与自变量之间的关系的方程叫做微分方 程．未知函数是一元函数的，叫常微分方 程;未知函数是多元函数的，叫做偏微分方 程．
如
y ' x y'x2y2 y''y'xy
Matlab实现 [t,x]=ode45(f,ts,x0,options,p1,p2,......)
50 0 0
30 20 10
0 0
10
20
50
30
20
10
0
30
0
10
8
6
4
2
100
0
50
100
150
50
100
高阶常微分方程的解法
高阶常微分方程
y ( n ) a 1 ( x ) y ( n 1 ) a ( n 1 ) ( x ) y ' a n ( x ) y f( x )
数值分析Ch7
2
f (xi , yi )
0
言
Numerical Analysis
M. H. Xu
二 理论基础 定理 7. 1 若f (x, y )在区 域 D = {(x, y )|a ≤ x ≤ b, |y | < ∞} 内连 续, 关于 y 满足 Lipschitz条件 , 即 存在 常数 L > 0, 对 任意 y1 , y2 , 不 等 式 |f (x, y1 ) − f (x, y2 )| ≤ L|y1 − y2 | 对所 有的 x ∈ [a, b]成立 , 则 初值 问题 dy = f (x, y ), a ≤ x ≤ b dx y (a) = y
Numerical Analysis
M. H. Xu
§7. 2 Euler方法与梯 形方 法 一 方法导出 由微分方程知 y (xi ) = f (xi , y (xi )), 用差商
y0 y yi y=y(x) yi+1 yn
y (xi+1 ) − y (xi ) h 近似 导数 y (xi )可 得 y (xi+1 ) ≈ y (xi ) + hf (xi , y (xi ))
0
在区 间[a, b]上有 唯 一解 y (x), 并且 y (x)为 连续 可微 的, 解函 数y (x) 连续 地依 赖于 初值 及f (x, y ).
Numerical Analysis
M. H. Xu
三 数值解法的基本步骤 第一 步: 把区 间[a, b]进行 划分 , 通 常进 行n等 分, 节点 xi = a + ih, i = 0, 1, 2 · · · , n, 其中 h = (b − a)/n; 第二 步: 求y (x)在节 点xi 处 函数 值y (xi )的近 似值 yi , 得 一列 表 函数 ; 第三 步: 根据 需要 可由 插值 方法求得 函数 y (x)在 x处的 近似 值, 或 由 列表 函数 求得 y = y (x)的近 似函 数. 说明 : 数值 解法 的关 键在 于如 何由 y0 得到y (x1 )的近 似值 y1 , 一 般地 , 如何 由y (xi )的近 似值 yi 得到y (xi+1 )的近 似值 yi+1 .
f (xi , yi )
0
言
Numerical Analysis
M. H. Xu
二 理论基础 定理 7. 1 若f (x, y )在区 域 D = {(x, y )|a ≤ x ≤ b, |y | < ∞} 内连 续, 关于 y 满足 Lipschitz条件 , 即 存在 常数 L > 0, 对 任意 y1 , y2 , 不 等 式 |f (x, y1 ) − f (x, y2 )| ≤ L|y1 − y2 | 对所 有的 x ∈ [a, b]成立 , 则 初值 问题 dy = f (x, y ), a ≤ x ≤ b dx y (a) = y
Numerical Analysis
M. H. Xu
§7. 2 Euler方法与梯 形方 法 一 方法导出 由微分方程知 y (xi ) = f (xi , y (xi )), 用差商
y0 y yi y=y(x) yi+1 yn
y (xi+1 ) − y (xi ) h 近似 导数 y (xi )可 得 y (xi+1 ) ≈ y (xi ) + hf (xi , y (xi ))
0
在区 间[a, b]上有 唯 一解 y (x), 并且 y (x)为 连续 可微 的, 解函 数y (x) 连续 地依 赖于 初值 及f (x, y ).
Numerical Analysis
M. H. Xu
三 数值解法的基本步骤 第一 步: 把区 间[a, b]进行 划分 , 通 常进 行n等 分, 节点 xi = a + ih, i = 0, 1, 2 · · · , n, 其中 h = (b − a)/n; 第二 步: 求y (x)在节 点xi 处 函数 值y (xi )的近 似值 yi , 得 一列 表 函数 ; 第三 步: 根据 需要 可由 插值 方法求得 函数 y (x)在 x处的 近似 值, 或 由 列表 函数 求得 y = y (x)的近 似函 数. 说明 : 数值 解法 的关 键在 于如 何由 y0 得到y (x1 )的近 似值 y1 , 一 般地 , 如何 由y (xi )的近 似值 yi 得到y (xi+1 )的近 似值 yi+1 .
常微分方程数值解法2线性多步法
对于线性多步法，其收敛性取决于微分方程的解的性质和方法的阶数。一般来说，高阶方法具有更好 的收敛性。
03
常见的线性多步法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最简单的一种方法，它基于线性近似，通过已知的函 数值来估计新的函数值。
详细描述
欧拉方法的基本思想是利用已知的函数值来估计下一个点的函数值。具体来说，假设我 们有一个函数 (y = f(x))，在已知 (x_0) 处的函数值 (y_0 = f(x_0)) 的情况下，欧拉方法 通过线性插值来估计 (x_1) 处的函数值 (y_1)，即 (y_1 = y_0 + h f(x_0))，其中 (h) 是
05
线性多步法的优缺点
优点
稳定性好
线性多步法在处理常微分方程时具有较好的数值稳定性， 能够有效地抑制数值振荡，提高计算结果的精度。
01
易于实现
线性多步法的计算过程相对简单，易于 编程实现，适合于大规模数值计算。
02
03
精度可调
通过选择不同的步长和线性多步法公 式，可以灵活地调整计算结果的精度， 满足不同的数值模拟需求。
改进方法的收敛性
研究收敛性条件
深入研究线性多步法的收敛性条件，了解哪些情况下方法可能不收 敛，并寻找改进措施。
优化迭代算法
通过优化迭代算法，提高方法的收敛速度和精度，减少迭代次数， 提高计算效率。
引入预处理技术
利用预处理技术对线性系统进行预处理，改善系统的条件数，提高方 法的收敛性。
拓展应用领域
在工程问题中的应用
控制系统设计
在工程领域中，线性多步法可以用于控制系统设计，通过 建立控制系统的数学模型，设计控制算法和控制器，实现 系统的稳定性和性能优化。
03
常见的线性多步法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最简单的一种方法，它基于线性近似，通过已知的函 数值来估计新的函数值。
详细描述
欧拉方法的基本思想是利用已知的函数值来估计下一个点的函数值。具体来说，假设我 们有一个函数 (y = f(x))，在已知 (x_0) 处的函数值 (y_0 = f(x_0)) 的情况下，欧拉方法 通过线性插值来估计 (x_1) 处的函数值 (y_1)，即 (y_1 = y_0 + h f(x_0))，其中 (h) 是
05
线性多步法的优缺点
优点
稳定性好
线性多步法在处理常微分方程时具有较好的数值稳定性， 能够有效地抑制数值振荡，提高计算结果的精度。
01
易于实现
线性多步法的计算过程相对简单，易于 编程实现，适合于大规模数值计算。
02
03
精度可调
通过选择不同的步长和线性多步法公 式，可以灵活地调整计算结果的精度， 满足不同的数值模拟需求。
改进方法的收敛性
研究收敛性条件
深入研究线性多步法的收敛性条件，了解哪些情况下方法可能不收 敛，并寻找改进措施。
优化迭代算法
通过优化迭代算法，提高方法的收敛速度和精度，减少迭代次数， 提高计算效率。
引入预处理技术
利用预处理技术对线性系统进行预处理，改善系统的条件数，提高方 法的收敛性。
拓展应用领域
在工程问题中的应用
控制系统设计
在工程领域中，线性多步法可以用于控制系统设计，通过 建立控制系统的数学模型，设计控制算法和控制器，实现 系统的稳定性和性能优化。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
251 5 ( 5 ) h y ( i ) 4阶Adams显式公式的截断误差为 720 19 5 ( 5 ) y ( x ) y = h y ( i ) 4阶Adams隐式公式的截断误差为 i +1 i +1 720 y( x ) - y 251 当 h 充分小时，可近似认为i i ，则： y( xi +1 ) - yi +1 - 19 i +1 i +1 y ( x i +1 ) - y i +1 =
y(xi+1). 其通式可写为：
yi +1 = a0 yi + a1 yi -1 + ... + ak yi -k + h(b-1 fi +1 + b0 fi + b1 fi -1 + ... + bk fi -k )
基于数值积分的构造法
将 y = f ( x , y ) 在 [ xi , xi +1 ] 上积分，得到
h yi +1 = yi + h [ f i + t ( f i - f i -1 )] dt = yi + ( 3 f i - f i -1 ) 0 2
1
Ri = h
1
0
d 2 f ( x , y( x )) 1 5 3 t h ( t + 1 ) h dt = h y ( i ) 2 dx 2! 12
0 0
1
1
yi +1 = yi + h N k (xi + t h) dt
0
1
/* 显式计算公式 */ Newton
1
局部截断误差为： Ri = y( xi +1 ) - yi +1 = h Rk ( xi + t h) dt
0
插值余项
例1 k=1 时有 N1 ( xi + t h) = f i + t f i = f i + t ( f i - f i -1 )
一、收敛性
定义 若某算法对于任意固定的 x = xi = x0 + i h，当 h0
( 同时 i ) 时有 yi y( xi )，则称该算法是收敛的. 例3
y = y 就初值问题 y (0) = y0
利用k+1 个节点上的被积函数值 fi+1 , fi , …, fi-k+1 构造 k 阶 牛顿前插多项式。与显式多项式完全类似地可得到一系列 ~ k + 2 ( k + 2) ~ 隐式公式，并有 Ri = Bk h y (i ) ，其中 Bk 与 fi+1 , fi , …, fi-k+1 的系数亦可查表得到.
5 12 37 24 -
Misprint on p.106
9 24
5 12
3 8
251 720
…
常用的是 k = 3 的4阶亚当姆斯显式公式
y i +1 h = yi + (55 f i - 59 f i -1 + 37 f i - 2 - 9 f i - 3 ) 24
上页 下页
返回
…
…
…式公式
y( xi +1 ) - y( xi ) =
xi +1 xi
f ( x, y( x))dx
xi +1
i
只要近似地算出右边的积分 I k f ( x, y( x ))dx ，则可通 x 过 yi +1 = yi + I k 近似y(xi+1) .而选用不同近似式 Ik，可得到不 同的计算公式.
251 ( y i +1 - y i +1 ) 270 19 y ( x i +1 ) y i +1 ( y i +1 - y i +1 ) 270 y ( x i +1 ) y i +1 +
外推技术 Modified final Corrected Modified /* extrapolation */ value valuey c i+1 value m i+1 i+1
上页 下页
返回
例2 给 出 一 个 三 步 显 式 格 式 y n+1 = y n + h(af n-1 + b f n- 2 ), 其 中f n = f ( x n , y n ). 试 用Taylor 展开原理确定 a、b ， 使 该 公 式 具 有 尽可能高的阶 .
上页 下页
返回
§5 收敛性和稳定性
上页 下页
返回
亚当姆斯显式公式
利用k+1 个节点上的被积函数值 f i , f i -1 , ... , f i -k 构造 k 阶牛顿 后插多项式 N k ( xi + t h) , t [0, 1], 有
xi +1 xi
f ( x, y( x))dx = N k ( xi + t h) h dt + Rk ( xi + t h) h dt
第七章 常微分方程的数值解法（二）
第四节 第五节 第六节 第七节 线性多步法 单步法的收敛性与稳定性 一阶方程组和高阶方程 边值问题的数值解法
§4 线性多步法
当 b-10 时，为隐式公 式; b-1=0 则为显式公式. 用若干节点处的 y 及 y’ 值的线性组合来近似 f j = f (xj , yj )
…
…
…
…
…
较同阶显 式稳定
上页 下页
返回
…
亚当姆斯预测-校正系统
Predicted 注意：三步所用公 value pi+1 Step 1: 用Runge-Kutta 法计算前 k 个初值； 式的精度必须相同。 通常用经典RungeStep 2: 用Adams 显式计算预测值； Kutta 法配合4阶 Adams 公式. Step 3: 用同阶Adams 隐式计算校正值.
k 0 1 fi+1 1
1 2 5 12 9 24 1 2 8 12 19 24
fi
f i- 1
f i- 2
…
Bk
1 2
~
小于Bk
1 12 5 24 1 24
-
1 12
2
3
1 24 19 720
…
常用的是 k = 3 的4阶亚当姆斯隐式公式
y i +1 = y i + h (9 f i +1 + 19 f i - 5 f i -1 + f i - 2 ) 24
上页 下页
返回
注：一般有 Ri = Bk hk +2 y( k +2) (i ) ，其中Bk 与yi+1 计算公式
中 fi , …, fi-k 各项的系数均可查表得到 .
k 0 fi 1
3 2 1 2
f i- 1
f i- 2
f i- 3
…
Bk
1 2
1
2 3
23 12 55 24
-
16 12 59 24
y(xi+1). 其通式可写为：
yi +1 = a0 yi + a1 yi -1 + ... + ak yi -k + h(b-1 fi +1 + b0 fi + b1 fi -1 + ... + bk fi -k )
基于数值积分的构造法
将 y = f ( x , y ) 在 [ xi , xi +1 ] 上积分，得到
h yi +1 = yi + h [ f i + t ( f i - f i -1 )] dt = yi + ( 3 f i - f i -1 ) 0 2
1
Ri = h
1
0
d 2 f ( x , y( x )) 1 5 3 t h ( t + 1 ) h dt = h y ( i ) 2 dx 2! 12
0 0
1
1
yi +1 = yi + h N k (xi + t h) dt
0
1
/* 显式计算公式 */ Newton
1
局部截断误差为： Ri = y( xi +1 ) - yi +1 = h Rk ( xi + t h) dt
0
插值余项
例1 k=1 时有 N1 ( xi + t h) = f i + t f i = f i + t ( f i - f i -1 )
一、收敛性
定义 若某算法对于任意固定的 x = xi = x0 + i h，当 h0
( 同时 i ) 时有 yi y( xi )，则称该算法是收敛的. 例3
y = y 就初值问题 y (0) = y0
利用k+1 个节点上的被积函数值 fi+1 , fi , …, fi-k+1 构造 k 阶 牛顿前插多项式。与显式多项式完全类似地可得到一系列 ~ k + 2 ( k + 2) ~ 隐式公式，并有 Ri = Bk h y (i ) ，其中 Bk 与 fi+1 , fi , …, fi-k+1 的系数亦可查表得到.
5 12 37 24 -
Misprint on p.106
9 24
5 12
3 8
251 720
…
常用的是 k = 3 的4阶亚当姆斯显式公式
y i +1 h = yi + (55 f i - 59 f i -1 + 37 f i - 2 - 9 f i - 3 ) 24
上页 下页
返回
…
…
…式公式
y( xi +1 ) - y( xi ) =
xi +1 xi
f ( x, y( x))dx
xi +1
i
只要近似地算出右边的积分 I k f ( x, y( x ))dx ，则可通 x 过 yi +1 = yi + I k 近似y(xi+1) .而选用不同近似式 Ik，可得到不 同的计算公式.
251 ( y i +1 - y i +1 ) 270 19 y ( x i +1 ) y i +1 ( y i +1 - y i +1 ) 270 y ( x i +1 ) y i +1 +
外推技术 Modified final Corrected Modified /* extrapolation */ value valuey c i+1 value m i+1 i+1
上页 下页
返回
例2 给 出 一 个 三 步 显 式 格 式 y n+1 = y n + h(af n-1 + b f n- 2 ), 其 中f n = f ( x n , y n ). 试 用Taylor 展开原理确定 a、b ， 使 该 公 式 具 有 尽可能高的阶 .
上页 下页
返回
§5 收敛性和稳定性
上页 下页
返回
亚当姆斯显式公式
利用k+1 个节点上的被积函数值 f i , f i -1 , ... , f i -k 构造 k 阶牛顿 后插多项式 N k ( xi + t h) , t [0, 1], 有
xi +1 xi
f ( x, y( x))dx = N k ( xi + t h) h dt + Rk ( xi + t h) h dt
第七章 常微分方程的数值解法（二）
第四节 第五节 第六节 第七节 线性多步法 单步法的收敛性与稳定性 一阶方程组和高阶方程 边值问题的数值解法
§4 线性多步法
当 b-10 时，为隐式公 式; b-1=0 则为显式公式. 用若干节点处的 y 及 y’ 值的线性组合来近似 f j = f (xj , yj )
…
…
…
…
…
较同阶显 式稳定
上页 下页
返回
…
亚当姆斯预测-校正系统
Predicted 注意：三步所用公 value pi+1 Step 1: 用Runge-Kutta 法计算前 k 个初值； 式的精度必须相同。 通常用经典RungeStep 2: 用Adams 显式计算预测值； Kutta 法配合4阶 Adams 公式. Step 3: 用同阶Adams 隐式计算校正值.
k 0 1 fi+1 1
1 2 5 12 9 24 1 2 8 12 19 24
fi
f i- 1
f i- 2
…
Bk
1 2
~
小于Bk
1 12 5 24 1 24
-
1 12
2
3
1 24 19 720
…
常用的是 k = 3 的4阶亚当姆斯隐式公式
y i +1 = y i + h (9 f i +1 + 19 f i - 5 f i -1 + f i - 2 ) 24
上页 下页
返回
注：一般有 Ri = Bk hk +2 y( k +2) (i ) ，其中Bk 与yi+1 计算公式
中 fi , …, fi-k 各项的系数均可查表得到 .
k 0 fi 1
3 2 1 2
f i- 1
f i- 2
f i- 3
…
Bk
1 2
1
2 3
23 12 55 24
-
16 12 59 24