数值分析21(常微分方程数值解)
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
12/27
基本概念
单步法: 即为了求得后一点xn+1上数值解yn+1,只要知 道一点xn的数值解yn就可以了。
多步法: 即计算yn+1时,要用到前面多点的数值解yn, yn-1, · · · , yn-r等信息。 显式格式:数值解yn+1可以用xn和yn 解析表出。这类 格式称为显式格式。反之称为隐式格式。对隐式格 式而言,计算数值解yn+1形式上需要解方程(组) 。 显式方法易于计算, 但通常数值稳定性较差。而隐 式方法通常稳定性较好, 但计算上显然很不方便。 将显式方法和隐式方法联合使用是重要的思想方法。
洛伦兹吸引子 lorenzgui
关于爱情的动力系统 Lover affaires and differential equations
常微分方程的初值问题是描述系统发展演变的重 要工具和手段。
2/27
Euler法与修正的Euler法
一阶常微分方程初值问题: dy f ( x , y ), x x 0 其中y = y(x) 是未知函数,
右矩阵公式:
xn1 xn
f ( x, y( x ))dx hf ( xn1 , yn1 )
隐式Euler法
yn1 yn hf ( xn1 , yn1 )
梯形公式:
x n 1 xn
h f ( x , y( x ))dx [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
dx y ( x 0 ) y0
右端函数 f(x, y )是已知函数, 初值 y0 是已知数据。 数值方法——取定离散点: x0 < x1 < x2 < · · ·< xN 求未知函数 y(x) 在离散点处的近似值 y1, y2, y3, · · · · · , yN
考察点xn列 出微分方程
dy 例1. 用Euler法求初值问题 y xy2 , 0 x 2 dx 的数值解。 y ( 0) 1
解: 记 f (x, y) = y－ x y2, xn= nh (n = 0, 1, 2,· · · ,N) 由Euler公式得: yn+1 = yn + h( yn－ xn yn2 ) (n = 0, 1, · · · ,N)
《数值分析》 23
解一阶常微分方程欧拉法
局部截断误差与p阶精度
Range-Kutta公式
一阶常微分方程组和二阶方程
线性多步法简介
1/27
常微分方程(ODE)的定解问题主要有初值问题 (IVP)和边值问题(BVP)两大类。
急性传染病数学模型(SIR模型)
di si i , i (0) i0 dt ds si , s(0) s0 dt
微积分基本定理 f ( x )dx F (b) F (a )
a
y( xn1 ) y( xn )
左矩形公式
xn
f ( x, y( x ))dx
xn1
xn
f ( x, y( x ))dx hf ( xn , yn )
显式Euler法
6/27
yn1 yn hf ( xn , yn )
2 1.5
10
0.2 0.0123 0.1059
20
0.1 0.0026 0.0521
1.5
30
0.0667 0.0011 0.0342
40
0.05 5.9612e-004 0.0256
1 1 0.5 0.5 0 -1 0 -1
0
1
2
3
4
5
0
1
2
3
4
5
欧拉方法(h=0.2) 误差最大值 0.1059
20
0.1 0.0026 0.0521
2 1.5
30
0.0667 0.0011 0.0342
40
0.05 5.9612e-004 0.0256
1 1 0.5 0.5 0 -1 0 -1
0
1
2
3
4
5
0
1
2
3
4
5
预-校方法, h=0.2时 误差最大值: 0.0123
欧拉方法,
h=0.2时
误差最大值: 0.1059
h yn1 yn [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
7/27
预报-校正方法(修正的Euler法):
h yn1 yn [ f ( xn , yn ) f ( x n1 , yn1 )] 2
~ yn1 yn hf ( xn , yn ) 预报
四阶Range-Kutta公式一般形式
yn+1= yn+ h[k1+2k2+2k3+k4]/6 k1=f(xn,yn), k2=f(xn+0.5h, yn+0.5hk1) k3=f(xn+0.5h, yn+0.5hk2), k4=f(xn+h, yn+hk3)
18/27
dy 2 y xy , 0 x 2 例4 dx 1 y( x ) x y ( 0 ) 1 x 1 2 e
4/27
取步长 h = 2/10, 2/20, 2/30, 2/40, 用Euler法求解 的数值实验结果如下.
N
h 误差
2 1.5 1 0.5 0 -1
10
0.2 0.1059
20
0.1 0.0521
30
0.0667 0.0342
40
0.05 0.0256
解析解:
0 1 2 3 4 5
o —— 数值解 ---- —— 准确解
i 1
s
取点个数s称为龙格-库塔法的级数。一般的RK公式 如下给出: K 1 f ( x n , yn ) K 2 f ( xn 2 h, yn h21 K 1 )
K s f ( xn s h, yn h( s1 K 1 yn1 yn h ci K i
预-校方法(h=0.2) 误差最大值0.0123
9/27
局部截断误差
Tn+1=y(xn+1) - yn+1称为局部截断误差。设 yn= y(xn)
由Taylor公式
仅考虑xn到xn+1的局部情况,并假定xn之前的计算没有误差
( xn1 xn ) 2 y( xn1 ) y( xn ) ( xn1 xn ) y( xn ) y( ) 2 2 h y( ) 即 y( xn1 ) y( xn ) hf ( xn , yn ) 2 Euler公式: yn+1 = yn+ hf (xn, yn) 的局部截断误差
数值实验:几种不同求数值解公式的误差比较
n h RK4 10 0.2 20 0.1 30 0.0667 40 0.05 2.186e-007
6.862e-005 3.747e-006 7.071e-007
RK3
RK2 Euler
0.0012
0.0123 0.1059
1.529e-004 4.517e-005 1.906e-005
y(xn+1) – yn+1=y(xn) – yn+ O(h2) = O(h2)
Euler公式的局部截断误差记为 O(h2) , 则称 Euler公式具有1阶精度。
10/27
从x0开始计算, 如果考虑每一步产生的误差直到xn,则 误差y(xn)-yn称为数值方法在节点xn处的整体误差。
粗略的讲, 整体误差把许多步的局部截断误差加在一起, 而步数正比于1/h(即步长倒数), 因此整体误差O(h p) 。
i 1 s
s , s 1 K s ))
RK方法是单步法(为了求得后一点xn+1上数值解yn+1, 只要知道一点xn的数值解yn)。
15/27
下面以2级方法为例子具体介绍龙格-库塔法。 yn1 yn hc1 K 1 hc2 K 2
K 1 f ( x n , yn ) K 2 f ( xn 2 h, yn h21 K 1 )
1 y( x ) x 1 2e x
5/27
求解常微分方程初值问题的数值方法可以由很多, 各 种方法都可以看成使用不同的数值积分方法计算 用数值积分方法离散化常微分方程 y' = f (x, y)
xn1
xn
b
y( x )dx
xn1
xn1
xn
f ( x, y( x ))dx
故方程组有无穷多组解, 每一组解构成的RK公式的阶数 都是2, 都叫做二阶RK公式。
c1 c2 1 / 2, 2 21 1就是改进Euler法。
17/27
三阶Range-Kutta公式一般形式
yn+1= yn+ h[k1+4k2+k3]/6 k1=f(xn,yn), k2=f(xn+0.5h, yn+0.5hk1) k3=f(xn+h, yn – hk1+2hk2)
yn1
h yn [ f ( xn , yn ) f ( x n1 , ~ yn1 )] 2
预-校方法算法如下 k 1 = f ( x n , yn ) , k2 = f( xn+1 , yn+ h k1) ,
校正
yn1
h y n [ k1 k 2 ] 2
8/27
n
h 误差2 Байду номын сангаас差1
构造的基本思想是选择适当的系数使得方法的局部 截断误差阶数尽可能高。
Tn1 y( xn1 ) yn hc1 f ( xn , yn ) hc2 f ( xn 2h, yn h21 f ( xn , yn ))
f ( xn 2 h, yn h21 f ( xn , yn )) f ( xn , yn ) 2 hf x ( xn , yn ) h21 f ( xn , yn ) f y ( xn , yn ) O( h2 )
y ( xn 1 ) y ( x n ) y( xn ) f ( xn , yn ) f ( x n , yn ) h
3/27
足够多的点来近 似连续对象
y( xn1 ) yn hf ( xn , yn ) 其中yn是y(xn)的近似。
求解常微分方程初值问题的Euler方法 取定步长: h,记 xn = x0 + nh, ( n = 1,2, · · · ,N) 称计算格式: yn+1 = yn + h f( xn, yn ) 为显式Euler公式。 对应的求初值问题数值解的方法称为Euler方法。
若局部截断误差为O(h p +1)则称显式单步法具有 p 阶精度 。
例 2. 证明梯形方法具有2阶精度 h yn1 yn [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
Tn1 y( xn1 ) yn1
hy( xn )
h2 2
h y( xn1 ) y( xn ) [ y( xn ) y( xn1 )] 2
16/27
Tn1 (1 c1 c2 ) f ( xn , yn )h (1 / 2 c2 2 ) f x ( xn , yn )h2 +(1 / 2 c2 21 ) f y ( xn , yn ) f ( xn , yn )h2 O( h3 )
c1 c2 1, c22 1 / 2, c2 21 1 / 2
13/27
改进的Euler方法这样理解: 它用xn和xn+1两个点 的斜率值K1和K2取算术平均作为平均斜率, 而xn+1处 的斜率值则利用已知信息xn来预报。
~ yn1 yn hf ( xn , yn )
h yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 2 是否可以推广改进的Euler方法？
h2 2
y( xn )
h3 3!
h y ( xn ) [ y( xn ) y( xn ) hy( xn ) 2
y( x n )] O( h4 )
h3 y( xn ) O( h4 ) 12
11/27
n
h 误差2 误差1
1.5
10
0.2 0.0123 0.1059
(Runge-Kutta)龙格-库塔法
RK方法是一大类的方法, 其基本思想是采用如下形式:
yn1 yn h ci f ( i , yi )
i 1
s
其中i [ xn , xn1 ]且1 2
s。
14/27
yn1 yn h ci f ( i , yi )
基本概念
单步法: 即为了求得后一点xn+1上数值解yn+1,只要知 道一点xn的数值解yn就可以了。
多步法: 即计算yn+1时,要用到前面多点的数值解yn, yn-1, · · · , yn-r等信息。 显式格式:数值解yn+1可以用xn和yn 解析表出。这类 格式称为显式格式。反之称为隐式格式。对隐式格 式而言,计算数值解yn+1形式上需要解方程(组) 。 显式方法易于计算, 但通常数值稳定性较差。而隐 式方法通常稳定性较好, 但计算上显然很不方便。 将显式方法和隐式方法联合使用是重要的思想方法。
洛伦兹吸引子 lorenzgui
关于爱情的动力系统 Lover affaires and differential equations
常微分方程的初值问题是描述系统发展演变的重 要工具和手段。
2/27
Euler法与修正的Euler法
一阶常微分方程初值问题: dy f ( x , y ), x x 0 其中y = y(x) 是未知函数,
右矩阵公式:
xn1 xn
f ( x, y( x ))dx hf ( xn1 , yn1 )
隐式Euler法
yn1 yn hf ( xn1 , yn1 )
梯形公式:
x n 1 xn
h f ( x , y( x ))dx [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
dx y ( x 0 ) y0
右端函数 f(x, y )是已知函数, 初值 y0 是已知数据。 数值方法——取定离散点: x0 < x1 < x2 < · · ·< xN 求未知函数 y(x) 在离散点处的近似值 y1, y2, y3, · · · · · , yN
考察点xn列 出微分方程
dy 例1. 用Euler法求初值问题 y xy2 , 0 x 2 dx 的数值解。 y ( 0) 1
解: 记 f (x, y) = y－ x y2, xn= nh (n = 0, 1, 2,· · · ,N) 由Euler公式得: yn+1 = yn + h( yn－ xn yn2 ) (n = 0, 1, · · · ,N)
《数值分析》 23
解一阶常微分方程欧拉法
局部截断误差与p阶精度
Range-Kutta公式
一阶常微分方程组和二阶方程
线性多步法简介
1/27
常微分方程(ODE)的定解问题主要有初值问题 (IVP)和边值问题(BVP)两大类。
急性传染病数学模型(SIR模型)
di si i , i (0) i0 dt ds si , s(0) s0 dt
微积分基本定理 f ( x )dx F (b) F (a )
a
y( xn1 ) y( xn )
左矩形公式
xn
f ( x, y( x ))dx
xn1
xn
f ( x, y( x ))dx hf ( xn , yn )
显式Euler法
6/27
yn1 yn hf ( xn , yn )
2 1.5
10
0.2 0.0123 0.1059
20
0.1 0.0026 0.0521
1.5
30
0.0667 0.0011 0.0342
40
0.05 5.9612e-004 0.0256
1 1 0.5 0.5 0 -1 0 -1
0
1
2
3
4
5
0
1
2
3
4
5
欧拉方法(h=0.2) 误差最大值 0.1059
20
0.1 0.0026 0.0521
2 1.5
30
0.0667 0.0011 0.0342
40
0.05 5.9612e-004 0.0256
1 1 0.5 0.5 0 -1 0 -1
0
1
2
3
4
5
0
1
2
3
4
5
预-校方法, h=0.2时 误差最大值: 0.0123
欧拉方法,
h=0.2时
误差最大值: 0.1059
h yn1 yn [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
7/27
预报-校正方法(修正的Euler法):
h yn1 yn [ f ( xn , yn ) f ( x n1 , yn1 )] 2
~ yn1 yn hf ( xn , yn ) 预报
四阶Range-Kutta公式一般形式
yn+1= yn+ h[k1+2k2+2k3+k4]/6 k1=f(xn,yn), k2=f(xn+0.5h, yn+0.5hk1) k3=f(xn+0.5h, yn+0.5hk2), k4=f(xn+h, yn+hk3)
18/27
dy 2 y xy , 0 x 2 例4 dx 1 y( x ) x y ( 0 ) 1 x 1 2 e
4/27
取步长 h = 2/10, 2/20, 2/30, 2/40, 用Euler法求解 的数值实验结果如下.
N
h 误差
2 1.5 1 0.5 0 -1
10
0.2 0.1059
20
0.1 0.0521
30
0.0667 0.0342
40
0.05 0.0256
解析解:
0 1 2 3 4 5
o —— 数值解 ---- —— 准确解
i 1
s
取点个数s称为龙格-库塔法的级数。一般的RK公式 如下给出: K 1 f ( x n , yn ) K 2 f ( xn 2 h, yn h21 K 1 )
K s f ( xn s h, yn h( s1 K 1 yn1 yn h ci K i
预-校方法(h=0.2) 误差最大值0.0123
9/27
局部截断误差
Tn+1=y(xn+1) - yn+1称为局部截断误差。设 yn= y(xn)
由Taylor公式
仅考虑xn到xn+1的局部情况,并假定xn之前的计算没有误差
( xn1 xn ) 2 y( xn1 ) y( xn ) ( xn1 xn ) y( xn ) y( ) 2 2 h y( ) 即 y( xn1 ) y( xn ) hf ( xn , yn ) 2 Euler公式: yn+1 = yn+ hf (xn, yn) 的局部截断误差
数值实验:几种不同求数值解公式的误差比较
n h RK4 10 0.2 20 0.1 30 0.0667 40 0.05 2.186e-007
6.862e-005 3.747e-006 7.071e-007
RK3
RK2 Euler
0.0012
0.0123 0.1059
1.529e-004 4.517e-005 1.906e-005
y(xn+1) – yn+1=y(xn) – yn+ O(h2) = O(h2)
Euler公式的局部截断误差记为 O(h2) , 则称 Euler公式具有1阶精度。
10/27
从x0开始计算, 如果考虑每一步产生的误差直到xn,则 误差y(xn)-yn称为数值方法在节点xn处的整体误差。
粗略的讲, 整体误差把许多步的局部截断误差加在一起, 而步数正比于1/h(即步长倒数), 因此整体误差O(h p) 。
i 1 s
s , s 1 K s ))
RK方法是单步法(为了求得后一点xn+1上数值解yn+1, 只要知道一点xn的数值解yn)。
15/27
下面以2级方法为例子具体介绍龙格-库塔法。 yn1 yn hc1 K 1 hc2 K 2
K 1 f ( x n , yn ) K 2 f ( xn 2 h, yn h21 K 1 )
1 y( x ) x 1 2e x
5/27
求解常微分方程初值问题的数值方法可以由很多, 各 种方法都可以看成使用不同的数值积分方法计算 用数值积分方法离散化常微分方程 y' = f (x, y)
xn1
xn
b
y( x )dx
xn1
xn1
xn
f ( x, y( x ))dx
故方程组有无穷多组解, 每一组解构成的RK公式的阶数 都是2, 都叫做二阶RK公式。
c1 c2 1 / 2, 2 21 1就是改进Euler法。
17/27
三阶Range-Kutta公式一般形式
yn+1= yn+ h[k1+4k2+k3]/6 k1=f(xn,yn), k2=f(xn+0.5h, yn+0.5hk1) k3=f(xn+h, yn – hk1+2hk2)
yn1
h yn [ f ( xn , yn ) f ( x n1 , ~ yn1 )] 2
预-校方法算法如下 k 1 = f ( x n , yn ) , k2 = f( xn+1 , yn+ h k1) ,
校正
yn1
h y n [ k1 k 2 ] 2
8/27
n
h 误差2 Байду номын сангаас差1
构造的基本思想是选择适当的系数使得方法的局部 截断误差阶数尽可能高。
Tn1 y( xn1 ) yn hc1 f ( xn , yn ) hc2 f ( xn 2h, yn h21 f ( xn , yn ))
f ( xn 2 h, yn h21 f ( xn , yn )) f ( xn , yn ) 2 hf x ( xn , yn ) h21 f ( xn , yn ) f y ( xn , yn ) O( h2 )
y ( xn 1 ) y ( x n ) y( xn ) f ( xn , yn ) f ( x n , yn ) h
3/27
足够多的点来近 似连续对象
y( xn1 ) yn hf ( xn , yn ) 其中yn是y(xn)的近似。
求解常微分方程初值问题的Euler方法 取定步长: h,记 xn = x0 + nh, ( n = 1,2, · · · ,N) 称计算格式: yn+1 = yn + h f( xn, yn ) 为显式Euler公式。 对应的求初值问题数值解的方法称为Euler方法。
若局部截断误差为O(h p +1)则称显式单步法具有 p 阶精度 。
例 2. 证明梯形方法具有2阶精度 h yn1 yn [ f ( xn , yn ) f ( xn1 , y( xn1 ))] 2
Tn1 y( xn1 ) yn1
hy( xn )
h2 2
h y( xn1 ) y( xn ) [ y( xn ) y( xn1 )] 2
16/27
Tn1 (1 c1 c2 ) f ( xn , yn )h (1 / 2 c2 2 ) f x ( xn , yn )h2 +(1 / 2 c2 21 ) f y ( xn , yn ) f ( xn , yn )h2 O( h3 )
c1 c2 1, c22 1 / 2, c2 21 1 / 2
13/27
改进的Euler方法这样理解: 它用xn和xn+1两个点 的斜率值K1和K2取算术平均作为平均斜率, 而xn+1处 的斜率值则利用已知信息xn来预报。
~ yn1 yn hf ( xn , yn )
h yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 2 是否可以推广改进的Euler方法？
h2 2
y( xn )
h3 3!
h y ( xn ) [ y( xn ) y( xn ) hy( xn ) 2
y( x n )] O( h4 )
h3 y( xn ) O( h4 ) 12
11/27
n
h 误差2 误差1
1.5
10
0.2 0.0123 0.1059
(Runge-Kutta)龙格-库塔法
RK方法是一大类的方法, 其基本思想是采用如下形式:
yn1 yn h ci f ( i , yi )
i 1
s
其中i [ xn , xn1 ]且1 2
s。
14/27
yn1 yn h ci f ( i , yi )