微分方程数值解问题复习题

1+ λh + 1 λ2h2 + 1 λ3h3 + 1 λ4h4 ≤ 1
2!
3!
4!
这就是格式的绝对稳定区域。
二．抛物型方程的差分方程
1.(1) 简述用差分方法求解抛物型方程初边值问题的数值解的一般步骤。
(2)
写出近似一阶偏导数
(
∂u ∂x
)
n m
的三种有限差分逼近及其误差阶,写出近似
(
∂2u ∂x2
（2）
(
∂u ∂x
)nm
的三种有限差分逼近为：
向前差商：
(
∂u ∂x
)nm
≈
un m+1
−
umn
h
向后差商：
(
∂u ∂x
)nm
≈
umn
− umn −1 h
中心差商：
(
∂u ∂x
)nm
≈
un m+1
−
un m−1
2h
假设 u 有所要求的各阶偏导数，由泰勒公式：
un m+1
−
umn
h
=
(
∂u ∂x
)
n m
解答提示：
常微分方程初值问题的近似解公式
∑ ⎧
⎪ ⎨
yn+1
=
yn
+
k +1 j=0
a*j Δ j
fn− j+1
⎪⎩ y0 , y1,⋅⋅⋅, yk
∫ ∫ 称为 Admas 内插公式，其中， a*j = (−1) j
0⎛
−1
⎜ ⎝
−t j
⎞ ⎟ ⎠
dt
=
1 j!
0
t(t +1) ⋅⋅⋅ (t +
−1
的数值解，给出 y(1) 的近似值，精确到小数点后三位。
解答提示：
f (x, y) = − y + x +1，问题的预报-校正格式为
⎧ ⎪⎪ ⎨
y(0) n+1
yn+1
⎪
= =
yn yn
+ +
hf (xn , yn ) = yn + 0.1(−xn + yn
1 2
h
⎣⎡
f
(xn ,
yn ) +
f
( xn+1 ,
+
1 λ4h4 4!
)
yn
+
(1 +
λh
+
1λ 2!
2h2
+
1 λ3h3 3!
+
1λ 4!
4h4
)δn
=
yn+1
+
(1 +
λh
+
1 λ2h2 2!
+
1 λ3h3 3!
+
1 4!
λ
4h4
)δ n
因此，
δn+1
=
(1 +
λh
+
1λ 2!
2h2
+
1 λ3h3 3!
+
1λ 4!
4h4
)δ n
为了使得格式绝对稳定，必须 δn+1 ≤ δn ，也就是
)
n m
的差分逼近及其误差阶。
解答提示：
(1) 先构造微分方程的有限差分逼近。为此，首先将求解区域 Ω 用二组平行于坐
5
标轴的网格覆盖，选定适当的网格边长，用差商近似地替代偏导数，从而得到相
应的差分方程。
考虑原问题的初值条件和边值条件，从而得到差分方程的初始条件，并得到
部分边界节点的函数值信息。
最后求解差分方程，其解就是微分方程的近似解。
)mn
h
+⋅⋅⋅
因此，向前差商的误差为 O(h) 。
un m+1
−
un m−1
=
⎡⎢umn ⎣
+
(
∂u ∂x
)mn
h
+
1 2
(
∂2u ∂x2
)mn
h
2
+ ⋅⋅⋅⎤⎥⎦
−
⎡⎢umn ⎣
−
(
∂u ∂x
)
n m
h
+
1 2
(
∂2u ∂x2
)mn
h2
+ ⋅⋅⋅⎤⎥⎦
2h
2h
=
( ∂u ∂x
)
n m
+
1 (∂3u 3! ∂x3
3.(1) 叙述 Admas 外插公式；
1
(2) 设 a j 为 Admas 外插公式的系数，计算 a4 的准确值。 解答提示：
常微分方程初值问题
的近似解公式
⎧⎪ dy ⎨ dx
=
f
(x,
y)
⎪⎩ y(x0 ) = y0
称为 Admas 外插公式。
∑ ⎧
⎪ yn+1
=
yn
+
k
ajΔ j fn− j
+
1 λ2h2 2!
+
1 λ3h3 3!
+
1 4!
λ
4
h
4
)
yn
给 yn 一个扰动 δn ，将 yn = yn + δn 代入，求得
y n +1
=
yn+1
+
δn+1
=
(1 +
λh
+
1 λ2h2 2!
+
1 λ3h3 3!
+
1 4!
λ
4
h
4
)(
yn
+
δn )
=
(1 +
λh
+
1 λ2h2 2!
+
1 λ3h3 3!
因此，
en+1
=
y(xn+1) −
y(xn )
−
h(c1K1*
+
c2
K
* 2
)
=
f
(xn ,
y(xn ))h +
1 2
⎡⎣
fx (xn ,
y(xn )) +
f y (xn ,
y(xn )) f
(xn , y(xn ))⎤⎦ h2
+ O(h3)
− c1hf (xn , y(xn )) − c2h ⎡⎣ f (xn , y(xn )) + fx (xn , y(xn ))a2h + f y (xn , y(xn )) f (xn , y(xn ))b21h + O(h2 )⎤⎦
一.常微分方程数值解问题 1.以 h = 0.1 为步长，用欧拉法求初值问题
⎧⎪ dy = xe− y − y ⎨ dx ⎪⎩ y(0) = 1
的数值解，给出 y(1) 的近似值，精确到小数点后三位。
解答提示：
f (x, y) = xe− y − y ，问题的欧拉格式为
⎧⎪ ⎨ ⎪⎩
yn+1 = y0 = 1
要求其局部截断误差为 O(h3 ) ，确定系数 c1, c2 , a2 , b21 所满足的条件，并具体地给 出一个二级二阶 Runge-Kutta 格式。 解答提示：
格式的局部截断误差为
en+1 = y(xn+1) − y(xn ) − h(c1K1* + c2K2* )
其中，
K1* ,
K
* 2
分别为
yn
+
hf
(xn ,
yn )
=
yn
+
0.1(xne− yn
−
yn )
按照这样的公式进行计算，直到得出 y10 ，将其作为 y(1) 的近似值。计算的步骤
是： y0 → y1 → y2 → ⋅⋅⋅ → y10 。 具体的计算过程请自己完成。
2.以 h = 0.1 为步长，用预报校正格式求解初值问题
⎧⎪ dy = − y + x +1 ⎨ dx ⎪⎩ y(0) = 1
j −1)dt 。
∫ a4*
=1 4!
0
t(t +1)(t + 2)(t + 3)dt
−1
2
具体的计算请自己进行。
5.二级二阶 Runge-Kutta 格式为
⎧ ⎪ ⎨
yn+1 = K1 = f
yn + (xn ,
h(c1K1 yn )
+
c2 K2
)
⎪⎩K2 = f (xn + a2h, yn + b21hK1)
h
+
1 2
(
∂2u ∂x2
)mn
h2
h
+⋅⋅⋅
=
(
∂u ∂x
)nm
+
1 2
(
∂2u ∂x2
)mn
h
+
⋅
⋅
⋅
因此，向前差商的误差为 O(h) 。
umn
−
un m−1
h
=
umn
−
⎣⎡⎢umn
−
(
∂u ∂x
)mn
h
+
h
1 2
(
∂2u ∂x2
)
n m
h2
+ ⋅⋅⋅⎤⎥ ⎦
=
(
∂u ∂x
)
n m
−
1 2
(
∂2u ∂x2
y(0) n+1
)⎦⎤
+ 1) = yn
+
0.05(− xn
+
yn
+1−
xn+1
+
y(0) n+1
+ 1)
⎪⎩ y0 = 1
按照这样的公式进行计算，直到得出 y10 ，将其作为 y(1) 的近似值。
计算的步骤是：
y0
→
y(0) 1
→
y1
→
y(0) 2
→
y2
→
⋅⋅⋅
→
y(0) 10
→
y10
具体的计算过程请自己完成。
1 −t(−t −1) ⋅⋅⋅ (−t − j +1)dt =
0
j!
1 j!
1
t(t +1) ⋅⋅⋅ (t +
0
j −1)dt
具体计算请自己进行。
∫ a4
=
1 4!
1
t(t +1)(t + 2)(t + 3)dt
0
4.(1) 叙述 Admas 内插公式；
(2) 设 a*j 为 Admas 内插公式的系数，计算 a4* 的准确值。
=
(1
−
c1
−
c2
)
f
(
xn
,
y(
xn
))h
+
(
1 2
−
a2c2
)
fx
(
xn
,
y ( xn
))h2
+
(
1 2
−
c2b21
)
f
y
( xn
,
y ( xn
))
f
(
xn
,
y ( xn
))h2
+
O(h3
)
3
⎧
⎧
根据
en+1
=
O(h3 )
，必须
⎪⎪⎪⎨⎪112−−ca1 2−cc2 2==00
，也就是
⎪⎪⎪⎨ca12c+2c=2 ⎪
考虑对试验方程 dy = λ y 运用这些格式。作为课程设计问题之一，具体的步 dx
骤已经在上课的时候讲过，请自己写上。例如，对于经典四级四阶 Runge-Kutta 格式，我们如此求其绝对稳定区域。
经典四级四阶 Runge-Kutta 格式为
⎧ ⎪
yn+1
⎪
=
yn
+
1 6
h(K1
+
2K2
+
2K3
+
1 6
h(K1
+
2K2
+
2K3
+
K4)
=
yn
+
1 6
h
⎡⎢⎣λ
yn
+
2(λ
+
1 2
λ 2h) yn
+
2(λ
+
1 2
λ2h
+
1 4
λ 3h2 )
yn
+
(λ
+
λ2h
+
1 2
λ 3h2
+
1 4
λ 4h3 ) yn
⎤ ⎥⎦
=
yn
+
1 6
h(6λ
+ 3λ 2h
+
λ3h2
+
1 4
λ 4h3 ) yn
=
(1 +
λh
+
⋅⋅⋅
因此，近似公式的误差为 O(h2 ) 。
⎛0 1
⎞
⎜ ⎜
1
0
1
⎟ ⎟
2.求 n 阶矩阵方阵 S = ⎜ % % % ⎟ 的特征值。
⎜ ⎜
1
0
1
⎟ ⎟
⎜⎝
1 0⎟⎠
解答提示：
此为对角矩阵，我们有如下结论：
⎛b c
⎞
⎜ ⎜
a
b
c
⎟ ⎟
n 阶三对角矩阵 A = ⎜ % % % ⎟ 的特征值为
K
* 2
=
f
( xn
+ a2h, y(xn ) + b21hK1*) =
f (xn , y(xn )) +
fx (xn , y(xn ))a2h +
f y (xn , y(xn ))b21hK1*
+
1[ 2
fxx (xn
+ θ a2h,
y(xn )
+ θ b21hK1*)a22h2
+
2
f xy
( xn
⎨
j=0
⎪⎩ y0 , y1,⋅⋅⋅, yk
∫ 其中， a j
= (−1) j
1⎛
0
⎜ ⎝
−t j
⎞⎟dt ⎠
就是
Admas
外插公式的系数。
⎛t ⎞
⎜ ⎝
j
⎟ ⎠
=
t(t
−1) ⋅⋅⋅ (t j!
−
j
+ 1)
，特别地，
⎛ ⎜ ⎝
t 0
⎞ ⎟ ⎠
=
1，
⎛t⎞ ⎜⎝1⎟⎠
=
t
。
∫ ∫ aj
= (−1) j
)mn
h
+
1 2!
(
∂2u ∂x2
)mn
h
2
+
1 3!
(
∂3u ∂x3
)mn
h3
+ ⋅⋅⋅⎤⎥⎦
−
2umn
+
⎡⎢⎣umn
−
(
∂u ∂x
)
n m
h
+
1 2!
(
∂2u ∂x2
)mn
h
2
−
1 3!
(
∂3u ∂x3
)mn
h3
+ ⋅⋅⋅⎤⎥⎦
h2
=
(
∂2u ∂x2
)
n m
+
2 ∂4u 4!( ∂x4
)mn
h2
=1 1。 2
⎪1 ⎪⎩ 2
−
c2b21
=
0
⎪⎩⎪b21c2
=
1 2
令 c1
=
c2
=
1 2
,
a2
=
b21
= 1，就得到了预报-校正格式：
⎧ ⎪
yn+1
⎪
=
yn
+
1 2
h(K1
+
K2 )
⎨K1 = f (xn , yn )
⎪ ⎪
K
2
=
f (xn
+ h, yn
+ hK1)
⎩
6.求二级二阶，三级三阶，四级四阶 Runge-Kutta 格式的绝对稳定区域。(分别选 一个作为代表，即可)(补考用) 解答提示：
+θ a2h,
y ( xn
)
+ θ b21hK1* )a2hb21hK1*
+ f y2 (xn + θ a2h, y(xn ) + θ b21hK1* )b221h2K1*2 ]
= f (xn , y(xn )) + fx (xn , y(xn ))a2h + f y (xn , y(xn )) f (xn , y(xn ))b21h + O(h2 )
K4 )
⎪K1 = f (xn , yn )
⎪⎪ ⎨
K
2
⎪
=
f
( xn
+
1 h, 2
yn
+
1 2
hK1 )
⎪ ⎪
K
3
⎪
=
f (xn
+
1 2
h,
yn
+
1 2
hK2
)
⎪⎩K4 = f (xn + h, yn + hK3 )
对于试验方程 dy = λ y ，其 f (x, y) = λ y ， dx
K1 = f (xn , yn ) = λ yn ，
K2
=
f (xn
+
1 2 h, yn
+
1 2
hK1
)
=
λ
(
yn
+
1 2
hK1
)
=
λ
yn
+
1 2
λ
hK1
= λ yn
+
1 2
λ
h
⋅
λ
yn
，
=
(λ
+
1 2
λ
2h)
yn
4
K3
=
f
( xn
+
1 2 h, yn
+
1 2
hK
2
)
=
λ
(
yn
+
1 2
hK
2
)
=
λ
yn
+
1 2
λ
hK2
= λ yn
+
1 λh(λ + 2
1 2
λ
2h)
yn
=
(λ
+
1 2
λ
2h
+
1 4
λ
3h2
)
yn
K4
=
f (xn
+ h, yn
+ hK3 ) = λ( yn
+ hK3) = λ yn
+ λhK3
= λ yn
+ λh(λ +
1 λ2h + 2
1 4
λ
3h2
)
yn
=
(λ
+
λ2h
+
1 2
λ3h2
+
1 4
λ
4h3 )
yn
因此，
yn+1
=
yn
+
)mn
h2
+
⋅⋅⋅
因此，中心差商的误差为 O(h2 ) 。
(
∂2u ∂x2
)
n m
≈
( ∂u ∂x
)n m+1
−
( ∂u ∂x
)
n m
h
≈
un m+1
−
umn
h
−
umn
−
un m−1
h
h
=
un m+1
−
2umn
+
un m−1
h2
6
un m+1
−
2umn
+
un m−1
h2
=
⎡⎢⎣umn
+
(
∂u ∂x
K1,
K2
中
yn
替换为
ቤተ መጻሕፍቲ ባይዱ
y(
xn
)
所得。
根据一元函数的泰勒公式以及 dy = f (x, y) ，有 dx
y(xn+1) −
y(xn )
=
f
(xn ,
y(xn ))h +
1 2
⎡⎣
fx (xn ,
y(xn )) +
f y (xn ,
y(xn )) f
(xn , y(xn ))⎤⎦ h2
+ O(h3)
K1* = f (xn , y(xn )) ，