CFD2020-第10讲-常微分方程数值解法

t x
y
z
x
y
z
FDM
空间离散
U Q(U)
FVM
时间推进 t
DG
推荐方法： Runge-Kutta法
2阶
U (1) U n tQ(U n )
U n1 1/ 2U n 1/ 2U (1) 1/ 2tQ(U (1) )
3 阶 （TVD型）
U (1) U n tQ(U n ) U (2) 3/ 4U n 1/ 4U (1) 1/ 4tQ(U ) (1) U n1 1/ 3U n 2 / 3U (2) 2 / 3tQ(U ) (2)
仅迭代一次(n=1),将
y(0) m 1
写入第二式
ym1
ym
h 2
f (xm , ym ) f (xm1, ym hf (xm , ym ))
改进的Euler方法
Copyright by Li Xinliang
7
P阶Taylor级数法
设初值问题的解具有p+1阶连续导数，利用Taylor公式
Euler公式
K1 f (xm , ym )
ym1 ym hf (xm , ym )
N 2
ym1 ym h b1K1 b2K2
K1 f (xm , ym )
K2 f (xm c2h, ym ha21K1)
• 将 K2 在(xm , ym ) 处展开
K2
f
h(c2
fx
a21K1 f y )
y'''(xm ) O(h4 )
y' f
利 用
y'' f x ff y
y''' fxx 2 ffxy f 2 f yy ( fx ffy ) f y
二阶：O(h3 )
b1 b2 1 b2c2 1/ 2
c2 a21
无穷多组解
b2a21 1/ 2
Copyright by Li Xinliang
ym1 ym h(xm , ym , h)
d dx
f
(x,
y(x))
fx
f fy
d2 dx2
f (x, y(x))
f xx 2 ff xy
f 2 f yy
fy( fx
ff y )
……
(*)式导数的计算比较复杂 离散情形（非解析形式）
（很少使用）
Copyright by Li Xinliang
向前Euler方法
ym1 ym hf (xm1, ym1)
向后Euler方法
ym1
ym
h 2
f (xm , ym ) f (xm1, ym1)
梯形法
迭求代解
y(0) m1
y(n) m1
ym ym
hf (xm , ym )
h 2
f (xm , ym )
f (xm1,
y ) (n1) m1
K1 f (xm , ym )
K2 f (xm c2h, ym ha21K1)
f
h(c2
fx
a21K1 f y )
h2 2
(c22
f xx
2c2a21K1 f xy
a221K12
fyy)
h3 6
(c23 f xxx 3c22a21K1 f xxy 3c2a221K12
f
xyy
y(xm1)
p j0
hj j!
y
(
j)
(
xm
)
(
h p
p1
1)!
y(
p1)
(m
)
y' f (x, y)
(x, y(x),h)
p h j1 j1 j!
d j1 dx j1
f (x, y(x))
(*)
y ( xm 1 )
y ( xm
)
h(xm ,
y(xm ),
h)
(
h p1 p 1)!
y( p1) (m )
f yyy)
m
O(h4 )
ym1 ym b1hK1 b2hK2 b3hK3 O(h4 )
ym b1hf b3h[ f
b2h[ h(c3 fx
f
h(c2 fx (a31K1 a32
a21K1 K2) fy
f )
y)
h2
2
h2 2
(c22
(c32 f xx
f xx 2c2a21K1 f xy 2c3 (a31K1 a32K2
• ym1 ym h(b1 b2 ) f h2b2 (c2 fx a21 ffy )
h3 2
b2 (c22
f xx
2c2a21
ff
xy
a221
f
2
f yy)
O(h4 )
• 将 y(xm1) 在 xm 处展开
y ( xm1 )
y(xm ) hy'(xm )
h2 2
y''(xm )
h3 3!
Rm1 y(xm1) y(xm ) h(xm , y(xm ), h)
从准确值y( xm )计算1步后得 到的值与准确值y(xm1) 之差
整体截断误差:
m1 y(xm1) ym1
P阶方法:
从初值出发，计算m+1步后 得到的值与准确值之差
局部截断误差 O(h p1) 显式方法:
右端项不包含 ym1 隐式方法:
xm1 xm
f (x, y(x))dx
• 利用Taylor公式
y ( xm1 )
y(xm ) hy'(xm )
h2 2
y''(m )
Copyright by Li Xinliang
5
单步法
计算ym1 的值只用到前一步的值 ym
ym1 ym h(xm , ym , h)
增量函数
局部截断误差:
h3 2
b2 (c22
f xx
2c2a21
ff
xy
a221
f
2
f yy)
O(h4 )
Copyright by Li Xinliang
10
Runge-Kutta方法
N阶
N
ym1 ym h bi Ki i 1
N 2 ym1 ym h b1K1 b2K2
K1 f (xm , ym ) K2 f (xm c2h, ym ha21K1)
4阶
U (1) U n 1/ 2tQ(U n ) U (2) U n 1/ 2tQ(U ) (1) U (3) U n tQ(U (2) ) U n1 1/ 3(U n U (1) 2U (2) U ) (3) 1/ 6tQ(U ) (3)
更高阶 ……
Copyright by Li Xinliang
右端项含有 ym1
• 一阶向后差商代替一阶导数
y(xm1) h
y(xm )
y' ( xm1 )
• 两边从xm 到xm1 积分,积分采用 右矩形公式
y(xm1) y(xm )
xm1 xm
f (x, y(x))dx
Copyright by Li Xinliang
6
改进的Euler方法
ym1 ym hf (xm , ym )
y
y' f (x, y)
y2
y1
y(x0 ) y0
初值条件
y0
因为f(x,y(x))恰好是待求精确解y(x)的斜率y’(x)
y(x)是在x-y平面上经过点(x0,y0)的曲线,该曲线
在点(x0,y0)处的切线的斜率为f(x0,y0).
y y(x)
x0
x1 x2
x
Euler方法
y y0 f (x0 , y0 )( x x0 )
x x1 x h
y1 y0 hf (x0 , y0 ) y2 y1 hf (x1, y1)
ym1 ym hf (xm , ym )
• 一阶向前差商代替一阶导数
y(xm1) h
y(xm )
y'(xm )
• 两边从 xm到xm1 积分,积分采用
左矩形公式
y(xm1) y(xm )
11
Runge-Kutta方法
N阶
N 2
ym1 ym h b1K1 b2K2
K1 f (xm , ym ) K2 f (xm c2h, ym ha21K1)
2阶格式满足的条件：
N
ym1 ym h bi Ki i 1
b1 b2 1 b2c2 1/ 2 b2a21 1/ 2
y Ce f (x)dx
y(x) c'(x) g(x)e f (x)dx
• 常微分方程的数值计算
P(x, y) Q(x, y) y' 0
U P, x
U Q y
d U (x, y(x)) 0 dx
U (x, y(x)) C
计算机的发展、绝大多数常微分方程难以求解析解
U f1(U) f2 (U) f3 (U) g1(U) g2 (U) g3 (U)
f (x, y) 1
(y x)
i 1
ci aij , i 2,, N j 1
新的ym1与原来的 x 在计 算过程中有相同的增量
N
bi 1
i 1
Copyright by Li Xinliang
13
N
Runge-Kutta方法
ym1 ym h bi Ki
i 1
N 3 将 K2 K3 在 (xm, ym ) 处展开
K1 f (xm , ym )
i 1
Ki f (xm cih, ym h aij K j ) j 1
bi , ci , aij 是待定常数
Copyright by Li Xinliang
9
Runge-Kutta方法
N阶
N
ym1 ym h bi Ki i 1
N 1 可定出: b1 1
无穷多组解
c2 a21
• b1 b2 1/ 2, a21 c2 1
（改进的Euler方法）
• b1 0, b2 1, a21 c2 1/ 2
（中点方法）
• b1 1/ 4, b2 3 / 4, a21 c2 2 / 3
（Heun方法）
Copyright by Li Xinliang
2
系统-----积分-微分方程组-----微分方程组------常微分方程
• 简单常微分方程的分析解
常系数线性常微分方程、可分离变量型方程、全微分方程等
分离变量
线性方程(常数变易）
全微分方程
y' f (x)g( y)
dy g( y)
f
(x)dx
C
y' f (x) y g(x)
y' f (x)y y C(x)e f (x)dx
……
Wm (t) U (m1) (t)
WW 12''
W2 W 3
Wm' f (t,W1,,Wm )
U (i) (t0 ) Ui0
i 0,, m 1
初值条件
Copyright by Li Xinliang
4
§ 10.1 单步法（Runge-Kutta方法）
假定常微分方程的初值问题是唯一可解的
a221K12 ) fxy
fyy) (a31K1
O(h3)] O(h4 ) a32K2 )2 f )
O(h3
)]
O(h4
)
a32[ f h(c2 fx a21K1 f y ) O(h2)]
a231K13
f
yyy
)
m
O(h4 )
K3 f (xm c3h, ym ha31K1 ha32K2 )
[f
h(c3
h3 6
(c33
fx
(a31K1
a32K2 )
fy)
h2 2
(c32
f xx
2c3 (a31K1
a32K2 )
f xy
(a31K1
a32 K2 )2
f yy)
f xxx 3c32 (a31K1 a32K2 ) f xxy 3c2 (a31K1 a32K2 )2 f xyy (a31K1 a32K2 )3
12
Runge-Kutta方法
N
ym1 ym h bi Ki i 1
Butcher表
K1 f (xm , ym )
i 1
Ki f (xm cih, ym h aij K j ) j 1
0
c2 a21 c3 a31 a32 cN aN1 aN ,N 1
b1 bN 1 bN
自动满足条件：
t x
y
z
x
y
z
U Q(U) t
Copyright by Li Xinliang
3
常微分方程（组）：
U ' f (t,U )
U (t0 ) U0
初值条件
m阶常微分方程:
U (m) f (t,U ',U '',..., U (m1) )
引入辅助函数：
W1(t) U (t) W2 (t) U (1) (t)
8
Runge-Kutta方法
基本思想：函数在一点的导数值可以用该点附近若 干点的函数值近似表示,因此可以把Taylor级数法中
的增量函数 改为 f 在这些点上函数值的组合,然后
用Taylor展开确定待定系数,使方法达到一定的阶。
N级Runge-Kutta方法:
N
ym1 ym h bi Ki i 1
计算流体力学讲义2020
第十讲 常微分方程数值解法
李新亮 申义庆
知识点： Runge-Kutta方法 线性多步法 刚性常微分方程的数值解法
1 Copyright by Li Xinliang
知识回顾：
NS方程的（空间）离散
U f1(U) f2 (U) f3 (U) g1(U) g2 (U) g3 (U)
h2 2
(c22
f xx
2c2a21K1 fxy
a221K12
f
yy
)
m
O(h3)
f
h(c2
fx
a21 ffy )
h2 2
(c22
f xx
2c2a21 f fxy
a221 f
2
f yy)
m
O(h3)
ym1 ym h(b1 b2 ) f h2b2 (c2 fx a21 ff y )