微分方程数值解第五章答案
第五章
⎧1, x
⎪∂u ∂u
+=0, u (x ,0) =⎨1/2, x =0, 试分别用左偏心格式、1. 对初值问题
∂t ∂x ⎪
x >0. ⎩0,
LW 格式计算其数值解u k , k =1,2,3,4, 取τ/h =1/2
.
解: 矩形网格剖分区域. 取空间步长h , 时间步长τ的矩形网格剖分区域, 用节点(j , k ) 表示坐标点(x j , t k ) =(jh , k τ) , j =0, ±1, ±2,...; k =0,1,2,3,4.
⎡∂u ⎤⎡∂u ⎤
（1）左偏心格式：⎢+⎢⎥=0，在t 上用向前差商，x 上用向后
⎣∂t ⎦j ⎣∂x ⎦j
差商，得
+1
−u k u k j j
k k
τ
+
k
u k j −u j −1
h
=0，因为τ/h =1/2，整理得到
+1
u k =j
1k 1
u j −1+u k j 22
⎧1, j
⎪0
把已知条件离散成u j =⎨1/2, j =0，则可以根据下一层求上一层的值得
⎪0, j >0⎩
到u k ，k =1，2，3，4，下图中节点处值即为求出来的u k 值：
中国地质大学(北京)廉海荣编 1
LW 格式： u
k +1
j
ar k a 2r 2k k k
=u −(u j +1−u j −1) +(u j +1+u k −2u j −1j ) 在
22
k j
本题中, a =1, r =τ/h =1/2, 整理得到：
+1u k j
⎧1, j
33k 1k ⎪0
u =，同理可根据边值条件=u k +u −u ⎨1/2, j =0，j j −1j j +1
848⎪0, j >0
⎩
根据下一层求上一层的值得到u k ，k =1，2，3，4，下图中节点处值即为求出来的u k 值：
中国地质大学(北京)廉海荣编 2
∂u ⎧∂u +a =0, 0
2. 试对初边值问题⎨u(x,0)=ϕ(x), 0≤x0建立以
⎪u(0,t)=ψ(t), 0≤t ≤T . ⎪⎩
下差分格式 （a ）（b ）
+1u k −u k j j
τ
+1u k −u k j j
+a
+1k +1
u k j +1−u j −1
2h
=0,
τ
试分析它们的稳定性。
k +1k +1k k
a u j +1−u j −1u j +1−u j −1+(+) =0. 22h 2h
解: 首先用矩形网格剖分区域. 取空间步长h , 时间步长τ的矩形网格剖分区域, 用节点(j , k ) 表示坐标点(x j , t k ) =(jh , k τ) , j =0,1, 2,...; ⎡T ⎤
k =0,1, 2, " ⎢⎥.
⎣τ⎦
(a) 对于初值问题的方程在(j , k ) 处取值, 然后用向后差分代替时间偏
中国地质大学(北京)廉海荣编
3
导，然后用中心差商代替空间偏导为，舍去误差项, 得到差分格式
+a =0， 2h τ
然后用k +1层代替k 层有差分格式
=0 (5.1) +a
2h τ
下面我们对这个差分格式的稳定性进行分析.格式等价于
1+1+1k +1k
u k +ar (u k j j +1−u j −1)=u j
2
k i σj h
令u k 代入上式，得 j =v e
+1
u k −u k j j
+1k +1
u k j +1−u j −1
k −1
u k j −u j
k
u k j +1−u j −1
1
v k +1e i σjh +ar (v k +1e i σjh e i σh −v k +1e i σjh e −i σh )=v k e i σjh
2
进一步简化得传播因子为：
G (σ, τ) =
1
1
1+ar (e i σh −e −i σh )2
=
1
1+iar sin σh
显然
G (σ, τ) =
11 ==
1+iar sin σh +iar sin σh 对于a >0，任意的r 可知G (σ, τ) 恒小于1，VonNeumann条件成立, 故差分格式无条件稳定。
(b) 对于初值问题的方程在(j , k ) 处取值, 然后再用向前差分代替时间偏导，接着用中心差商代替空间偏导数, 舍去误差项定义差分格式
+1u k −u k j j
k
u k j +1−u j −1
τ
+a
2h
=0， (5.2)
1
然后(5.1)和(5.2)式两边均乘以，做和有
2
中国地质大学(北京)廉海荣编
4
第五章 双曲型方程的数值解
+1u k −u k j j
- 5 -
τ
k +1k +1k k
a ⎛u j +1−u j −1u j +1−u j −1⎞
=0. +⎜+⎟⎟2⎜2h 2h ⎝⎠
下面对这个差分格式的稳定性进行分析,差分格式等价于
11+1+1k +1k k
u k +ar (u k ar (u k j j +1−u j −1)=u j −j +1−u j −1)
44
k i σj h
令u k 代入上式，得 j =v e
11
v k +1e i σjh +ar (v k +1e i σjh e i σh −v k +1e i σjh e −i σh )=v k e i σjh +ar (v k e i σjh e i σh −v k e i σjh e −i σh )44
继续化简有
⎡1⎤k +1⎡1⎤k
⎢1+2iar sin σh ⎥v =⎢1−2iar sin σh ⎥v ⎣⎦⎣⎦
1
1−iar sin σh 从而得到传播因子为G (σ, τ) =.
1
1+iar sin σh 2
显然
11
1−iar sin σh −iar sin σh
2G (σ, τ) ==1 ==
111+iar sin σh 1+iar sin σh 22VonNeumann 条件成立, 故差分格式无条件稳定。
∂u ⎧∂u
⎪∂t +a (x ) ∂x =0, 0
的边初值问题建3. 对变系数方程⎨u (x ,0) =ϕ(x ), 0
⎪u (0,t ) =ψ(t ), 0≤t ≤T ⎪⎩
立左偏心格式及LF 格式，并研究它们的稳定性。
中国地质大学(北京)廉海荣编
5
解：首先用矩形网格剖分区域. 取空间步长h , 时间步长τ的矩形网格剖分区域, 用节点(j , k ) 表示坐标点(x j , t k ) =(jh , k τ) , j =0,1, 2,...; ⎡T ⎤
k =0,1, 2, " ⎢⎥.
⎣τ⎦
下面建立左偏心格式, 方程在(j , k ) 处取值, 关于时间和空间的一阶偏
导数均用向前差商代替, 舍去误差项可定义差分格式
+1u k −u k j j
τ
令r =
+a j
k
u k j −u j −1
h
=0
τ
h
，以上等式可化为
+1k
u k =(1−a j r ) u k j j +a j r u j −1
⎧u 0=ϕj , j =0,1, 2, " ⎪
对于初边值问题可定义差分方程组为⎨k j
==" u ψ, k 0,1, 2, ⎪k ⎩0
k i σj h
令u k 代入左偏心格式, 有 j =v e
e i σjh v k +1=(1−a j r) e i σjh v k +a j r e i σjh e −ijh v k
由此可得传播因子为G (σ, τ) =(1−a j r) +a j r e −ijh . 显然
G (σ, τ) =(1−a j r +a j r cos σh ) 2+(a j r sin σh ) 2 =1+[2(a j r ) −2a j r ](1-cos σh )
2
2
2
即(a j r ) 2≤a j r 时G (σ, τ) ≤1. 又因为0
rC 1≤1.
下面建立LF 格式, 首先采用特征线法, 如图所示
中国地质大学(北京)廉海荣编
6
t =t k +1
D BC
t =t k
为了求方程的解在P 点的函数值, 利用特征线法可知u (P ) =u (D ) , 故只需求解函数在D 点的函数值, ADC三点时间坐标相等, 故可以把此三点方程
的解看成空间的一元函数, 可以利用AC 两点函数值的线性插值来求D 点的函数值, 记
P :(j , k +1) ⇒(x j , t k +1) =(jh ,(k +1) τ) A :(j −1, k ) ⇒(x j −1, t k ) =((j −1) h , k τ) B :(j , k ) ⇒(x j , t k ) =(jh , k τ) C :(j +1, k ) ⇒(x j +1, t k ) =((j +1) h , k τ)
对A 、C 点线性插值有u (x D , t k ) =
+1u k =j
x C −x D x −x A
u (x A , t k ) +D u (x C , t k ) , 即 x A −x D x C −x A
h −a j r 2h
u k j +1+
h +a j r 2h
1k
[(1−a j r ) u k u k j −1=j +1+(1+a j r ) u j −1]. 2
也可以用差商代替偏导数建立LF 格式. 关于时间的一阶偏导数用向前差商代替, 关于空间的一阶偏导数用中心差商代替, 舍去误差项可定义中心差分格式
+1u k −u k j j
τ
+a j
k
u k j +1−u j −1
2h
=0, 此等价于
r +1k
u k a j (u k =u k j j +j +1−u j −1)
2
其中r =
τ
h
, 用u k j =
1k
(u j +1+u k j −1) 代替代入上式得LF 格式. 2
7
中国地质大学(北京)廉海荣编
k i σj h
令u k 代入LF 格式, 可得传播因子为 j =v e
G (σ, τ) =a j r cos σh −ia j r sin σh
22易得G (σ, τ) =1−(1−a 2j r )sin σh , 容易看出a j r ≤1时G (σ, τ) ≤1. 又因
2
2
为0
4. 研究方程
∂u ∂u
+a =0, −∞
r a =1时的稳定性，这里r =τ/h 。
解：首先构造leap-frog 格式, 用矩形网格剖分区域. 取空间步长h , 时间步长τ的矩形网格剖分区域, 用节点(j , k ) 表示坐标点(x j , t k ) =(jh , k τ) , ⎡T ⎤
j =0, ±1, ±2,...; k =0,1, 2, " ⎢⎥. 方程在(j , k ) 处取值, 分别用一阶中心差
⎣τ⎦
商代替一阶偏导数导数，则可以定义差分格式有：
+1−1−u k u k j j
2τ
此等价于u j
k +1
+a
k
u k j +1−u j −1
2h
=0，
−1k
=u k −ar (u k j j +1−u j −1) ，r =τ/h .
该方程为三层格式，将其转换为双层格式组讨论其稳定性. 令
k −1
v k j =u j ，那么得到方程组：
中国地质大学(北京)廉海荣编 8
k +1k k k
⎧⎪u j =v j −ar (u j +1−u j −1)
， ⎨k +1k
⎪⎩v j =u j
⎛u k ⎞j
又令w =⎜k ⎟=q k e i σjh 代入上式得则得到传播矩阵为：
⎜v ⎟⎝j ⎠
k j
⎛−are i σh +are −i σh
G (σ, τ) =⎜
1⎝1⎞⎛−2ari sin σh 1⎞
⎟=⎜⎟ 10⎠0⎠⎝
传播矩阵
G (σ, τ) 的特征方程λ2+2ari sin σh λ−1=0有根
λ1,2
==−ari sin σh ±
当1−a 2r 2sin 2σh ≥0，即ar ≤1时，得到共轭的两个特征值λ1、λ2，且有λ1=λ2=a 2r 2sin 2σh +1−a 2r 2sin 2σh =1. 传播矩阵的谱半径
ρ(G (σ, τ)) =max λj =1，
j
2
2
故V on Neumann条件成立。下面分析矩阵族G (σ, τ) 是否一致有界.
事实上当r a =1时，差分格式是不稳的，那么不妨设a >0，则有
{
k
}
ra =1，此时得到相应的传播矩阵为G (σ, τ) =⎜
σh =
π
2
−2i +n π时有G (σ, τ) =⎛⎜
⎝1
⎛−2i sin σh 1⎞
⎟, 取10⎝⎠
1⎞，下面用归纳法求G k (σ, τ) .
⎟0⎠
中国地质大学(北京)廉海荣编 9
⎛−2i 1⎞⎛−2i 1⎞⎛−3
G 2(σ, τ) =⎜⎟⎜⎟=⎜
1010⎝⎠⎝⎠⎝−2i ⎛−3−2i ⎞⎛−2i 1⎞⎛4i
G 3(σ, τ) =⎜⎟⎜⎟=⎜
−2i 110⎝⎠⎝⎠⎝−3⎛4i −3⎞⎛−2i 1⎞⎛5
G 4(σ, τ) =⎜⎟⎜⎟=⎜
−3−2i 10⎝⎠⎝⎠⎝4i "
⎛−(k +1) −ki ⎞
G k (σ, τ) =i k ⎜⎟
k −1⎠⎝−ki
−2i ⎞
⎟1⎠−3⎞⎟−2i ⎠4i ⎞⎟−3⎠
k =1, 2,3"
显然G k (σ, τ) 关于k 无界，那么矩阵族{G k (σ, τ) }不是一致有界，则当
ra =1差分格式不稳定。 当a
⎛2i sin σh 1⎞π−2i 1⎞，σ=−+n π，同样得到G (σ, τ) =⎛h ，此时取G (σ, τ) =⎜⎟⎜⎟10⎠2⎝⎝10⎠
那么同理说明此时差分格式仍不稳定。
综上所述，构建的跳蛙（leap-frog ）格式在r a =1时不稳定。 下面举一个实际的例子来分析该格式的稳定性情况： 设系数
a =1， 并设初值和边值均为常数函数，不妨假定二者分别为0.1和0，
设x 及t 的范围都为0~1，那么相应的得到方程为：
⎧∂u ∂u
=0, ⎪+
⎨∂t ∂x ⎪⎩ϕ(x , 0) =0. 1,
0≤x ≤1, 0≤t ≤1
φ(0, t ) =0
首先需要说明以下两个问题：
（1） 由于三层格式需要预先知道第0层和第1层的值，而这里只给定
第0层的初值，故利用左偏心格式求出第1层的值。已经证明左
偏心格式在ar ≤1, (a >0) 的情况下是稳定的。
中国地质大学(北京)廉海荣编
10
（2） 对每一层的最后一个点，即N+1号点亦利用左偏心格式求解。
编程求得不同网格比r 下得到的图形：
分析：由于设边值为0，又a>0，那么根据特征性的性质随着t 的增加函数u 的值应当趋近于0，显然在r=1/4,1/2都是满足这一性质的。当r=2时，u 的值在后期发生严 畸变，这说明当r>1时差分格式并不稳定。
实际上这里要着重讨论的是r=1时的情况，但很可惜在所举例子中这时本应当表现出的不稳定性并不是十分突出，甚至可以说并没有体现。经过分析，认为这可能是由于初、边值的取值不当造成的，另外由于第一层及每层最后一个点的计算采用了左偏心格式，这也是可能导致错误的一个
中国地质大学(北京)廉海荣编
11
地方。当然，也不排除程序的编写错误。
附程序： clear;
h=0.01;tao=0.03; fain=0.1;edge=0; a=1;
xmin=0;xmax=1;tmin=0;tmax=1;
r=tao/h;N=fix(xmax/h);T=(tmax/tao);
r0=1/4; %这里特别设置r0，使满足左偏心格式稳定的条件 value(1,1:N+1)=fain; %赋初值 value(1:T+1,1)=edge; %赋边值 for j=2:N+1
value(2,j)=r0*a*value(1,j-1)+(1-r0*a)*value(1,j); %利用左偏心格式设置第一层的数值，由于编程的方便，这里下标设置为2 end
for k=2:T for j=2:N
value(k+1,j)=value(k-1,j)-a*r*(value(k,j+1)-value(k,j-1)); %利用leap-frog格式计算 end
value(k+1,N+1)=r0*a*value(k,N)+(1-r0*a)*value(k,N+1); %利用左偏心格式计算每层第N+1个点 end
[X,Y]=meshgrid(xmin:h:xmax,tmin:tao:tmax); mesh(X,Y,value); xlabel('x'); ylabel('t');
zlabel('u(x,t)'); hold on;
中国地质大学(北京)廉海荣编
12
surf(X,Y,value)
2
∂2u 2∂u =0, a 为正常数的加权格式的稳定性. 5. 试研究波动方程2−a 2∂t ∂x
1
2
k a 22τ2
δt
u j
=h 2[θδx
u k +1+(1−2θ) δ2k 2k −1j x u j +θδx u j ] 解: 将二阶中心差分展开, 有
1
k −1
k k +1a 2+1k +1+1
τ
2
(u j
−2u j
+u j
) =h 2[θ(u k j −1−2u j +u k j +1) +(1−2θ)(u k u k k k −1k −1−1
−1−2j +u j +1) +θ(u j −1−2u j +u k j j +1)]令v k +1
=u k j j , 那么格式等价于
⎧u k +1k 22k +1k +1+1
⎪−2u k v r [θ(u j −1−2u j +u k j j +j =a j +1) ⎪⎨ +(1−2θ)(u k k k ) +θ(v k k k
j −1−2u j +u j +1j −1−2v j +v j +1)] ⎪⎪⎩v k +1k j =u j
令w k k k ) Τ
j =(u j , v j ，得
2w
k +1
+⎛⎜−a 2r 2θ0⎞j
⎟w k +1⎛2a 2r θ0⎞k +1
⎝00⎠j +1+⎜⎝00⎟⎠
w j +⎛⎜
−a 2r 2θ0⎞k +1⎛a 2r 2⎝00⎟
⎠w (1−2θ) a 2r 2θ⎞k j −1=⎜⎝00⎟⎠
w j +1 +⎛⎜2−2a 2r 2(1−2θ) −1−2a 2r 2θ⎞k ⎛a 2r 2(1−2θ) a 2r 2θ⎞k
⎝10⎟⎠w j +⎜⎝00⎟⎠
w j −1
令 w k j
=q k e i σjh , 那么有q k +1=G (σ, τ) q k , 其中传播因子 中国地质大学(北京)廉海荣编
13
⎛2−2a 2r 2(1−2θ) +2a 2r 2(1−2θ) cos σh −1−2a 2r 2θ+2a 2r 2θcos σh ⎞
G (σ，τ) =⎜⎟
10⎝⎠1⎛⎞
0⎟⎜2222
⋅1+2a r θ−2a r θcos σh ⎜⎟⎜⎟01⎝⎠
⎛2−2a 2r 2(1−2θ) +2a 2r 2(1−2θ)cos σh ⎜
1+2a 2r 2θ−2a 2r 2θcos σh =⎜
1⎜
⎜
1+2a 2r 2θ−2a 2r 2θcos σh ⎝
⎞
−1−2a 2r 2θ+2a 2r 2θcos σh ⎟
⎟⎟
0⎟
⎠
其特征值方程为
2
2−2a 2r 2(1−2θ) +2a 2r 2(1−2θ)cos σh
λΙ−G (στ) =λ−λ+1=0
1+2a 2r 2θ−2a 2r 2θcos σh
特征值按模不大于1的充要条件是
2−2a 2r 2(1−2θ) +2a 2r 2(1−2θ)cos σh b =≤1−c =2
1+2a 2r 2θ−2a 2r 2θcos σh
成立, 而此充分条件为
⎧2−2a 2r 2(1−2θ) +2a 2r 2(1−2θ)cos σh ≤2+4a 2r 2θ−4a 2r 2θcos σh
⎨22222222
θθσθθσ−−+−≤−−+22a r (12) 2a r (12)cos h 24a r 4a r cos h ⎩
第一个不等式等价于−2a 2r 2+2a 2r 2θcos σh ≤0, 故对于任意的参数其都成立, 而第二式化简为2a 2r 2(1−cos σh )(4θ−1) +4≥0, 当θ≥等式成立, 而当θ∈(0,) 时, 若a 2r 2≤
1
4
1
时, 此不4
2
即ar ≤
(1−
cos σh )(1−4θ)
时, 不等式成立. 进一步可以证明此就是格式稳定的充分条件.
中国地质大学(北京)廉海荣编 14
∂w ⎧∂v a =⎪⎪∂t ∂x
6. 对方程组⎨，试建立LF 格式，并研究格式的稳定性.
w v ∂∂⎪=a ⎪∂t ∂x ⎩
解：用矩形网格剖分区域. 取空间步长h , 时间步长τ的矩形网格剖分区域, 用节点(j , k ) 表示坐标点(x j , t k ) =(jh , k τ) , j =0,1, 2,...; k =0,1, 2, " . 方程在(j , k ) 处取值, 关于时间的一阶偏导数用向前差商代替, 关于空间的一阶偏导数用中心差商代替, 舍去误差项可定义中心差分格式
+1k
⎧v k −v k w k j j j +1−w j −1
=a j ⎪
⎪τ2h
⎨k +1k k k
v j +1−v j −1⎪w j −w j
=a j ⎪τ2h ⎩
此等价于
ar k ⎧k +1k
=+(w j +1−w k v v j j j −1) ⎪⎪2 ⎨
ar ⎪w k +1=w k +(v k −v k )
j j j +1j −1⎪⎩2
其中r =
τ
h
, 用v k j =
1k 1k k k
(v j +1+v k ), (w =w +w j −1j j +1j −1) 代替代入上式22
得LF 格式.
ar k ⎧k +11k k
=++v (v v ) (w j +1−w k j j +1j −1j −1) ⎪⎪22 (1) ⎨
1ar ⎪w k +1=(w k +w k ) +(v k −v k )
j −1j j +1j −1j +1
⎪⎩22
⎛v ⎞τ
设m k j =⎜⎟, r =，(1)变为
h ⎝w ⎠
中国地质大学(北京)廉海荣编
k
15
第五章 双曲型方程的数值解
+1m k =j
- 16 -
ar ⎛01⎞k 1⎛10⎞k k k
⎜⎟(m j +1+m j −1) +⎜⎟(m j +1−m j −1) 2⎝01⎠2⎝10⎠
k i σj h
令m k ，则有 j =q e
q k +1e i σjh =
1⎛10⎞i σh ar ⎛01⎞i σh −i σh k i σjh −i σh k i σjh
() ++e e q e ⎜⎟⎜⎟(e −e ) q e
2⎝01⎠2⎝10⎠
iar sin σh ⎞k i σjh ⎛cos σh
=⎜⎟q e
cos σh ⎠⎝iar sin σh
故其LF 格式的传播因子为
iar sin σh ⎞⎛cos σh
G L (σ, τ) =⎜⎟， iar sin σh cos σh ⎝⎠
于是G L (σ, τ) =cos 2σh +a 2r 2sin 2σh =1−(1−a 2r 2)sin 2σh , 所以，当ar ≤1时，G L (σ, τ) ≤1，LF 格式稳定。
7. 试研究二维波动方程的显格式的稳定条件
k k −1u i k j +1=2(1−2r 2)u ij +r 2u i k +1, j +u i k −1, j +u i k , j +1+u i k , j −1−u ij 。
()
解:这是一个二维三层格式，利用传播因子法分析格式的稳定性。首先将其化为等价的双层格式组，有 u
k +1
(x , y )=2(1−2r 2)u k (x , y )+ 2
中国地质大学(北京)廉海荣编 16
第五章 双曲型方程的数值解 - 17 -
2⎧∂2u ∂u 2∂u =a −2C , 0
∂t ∂∂t x ⎪
∂⎪
8. 考虑带阻尼的波动问题⎨u (x ,0)=ϕ(x ), u (x ,0)=ψ(x ), 0≤x ≤1
∂t ⎪
⎪u (0, t )=u (1, t )=0, 0≤t ≤T ⎪⎩
其中a >0, C >0为常数。试证对上述方程建立的差分格式
δ2u k 222u k k +1
t j =a r δx j −2τC (u j −u k j )
当ar
解: 方法一: 类似于教材P200页讨论,
方法二: 类似于抛物型多层差分格式稳定性讨论 方法一: 1 把题中方程组化成一阶方程组 令υ=
∂u
, ω=a ∂u
∂t ∂x
， 则 ⎧⎪∂υ∂⎪⎨
∂t =a ω∂x
−2C υ ⎪∂ω∂υ ⎪⎩
∂t =a ∂x 对(5.3)建立显格式
⎧⎪υk +1−υk j ⎪=a ωk k j j +1/2−ωj −1/2
−2C ⎨τh
υk j
⎪ωk +1k k ⎪j −1/2−ωj −1/2υ+1j −υk +1
j −1
⎩
τ=a h 中国地质大学(北京)廉海荣编 (5.3) (5.4) 17
第五章 双曲型方程的数值解
k −1
u k j −u j
- 18 -
k
u k j −u j −1
若令υ=
k
j
τ
, ω
k
j −1/2
=a
h
, 由此可以推知格式(5.4)等价于题
目中的显格式，因此显格式的稳定性讨论可以通过格式(5.4)来进行。 把(5.4)整理变形, 有，
+1k k k k ⎧υk
⎪j −υj =ar (ωj +1/2−ωj −1/2)−2C τυj
⎨k +1
k k +1k +1−=−ωωar υυ (j j −1)⎪j −1/2⎩j −1/2
即
k +1k k k
⎧⎪υj =ar ωj +1/2+(1−2C τ) υj −ar ωj −1/2
(5.5) ⎨k +1k +1k +1k
ar υωar υω −++=⎪j j −1/2j −1j −1/2⎩
⎛υk ⎞j
令 M =⎜⎜ωk ⎟⎟, 那么(5.3)写成矩阵形式为 j ⎝⎠
k
j
⎛1⎜⎝−ar ⎛0=⎜⎝00⎞k +1⎛00⎞k +1⎛00⎞k +1⎟M j +⎜⎟M j −1/2+⎜⎟M j −1
ar 0⎠010⎝⎠⎝⎠
(5.6)
ar ⎞k ⎛1−2C τ0⎞k ⎛0−ar ⎞k
⎟M j +1/2+⎜⎟M j +⎜⎟M j −1/20⎠0001⎝⎠⎝⎠
k k i σjh
令M j =q e , 那么(5.6)式为
1⎛⎜
⎜−ar +are −i σh ⎝
0⎞k +1⎛1−2C τ⎟q =⎜−i σh ⎟⎜0e ⎠⎝
are
i σ−are
−i σ−i σe
⎞
⎟q k ⎟⎠
中国地质大学(北京)廉海荣编 18
q k +1
1⎛
=⎜
⎜−ar +are −i σh ⎝
−10⎞⎛⎜1−2C τi σ⎟⎜e ⎟⎠⎜0
⎝
ar e
(
i σ−e
−i σh
e
−i σ)⎞⎟q
⎟⎟⎠
k
iA ⎞k ⎛1−2C τ
=⎜q 2⎟C iA A (12τ) 1−−⎝⎠
其中A =2ar sin σ, 所以传播矩阵为
iA ⎞⎛1−2C τ
G (σ, r )=⎜2⎟ τ1−2C iA 1−A ()⎝⎠
4 而这个矩阵的特征方程为λ2−(2−2C τ−A 2)λ+1−2C τ=0 它的特征根为
λ1,2=
22
2
=1−C τ−2a r sin σ() (1) 若1−C τ−2a 2r 2sin 2σ−(1−2C τ) ≤0, 那么特征值为复根,
|λ1,2|2=1−C τ−2a 2r 2sin 2σ()
2
(
)
2
+(1−2C τ) −1−C τ−2a 2r 2sin 2σ(
)=1−2C τ
满足von Neumann条件. 此时显然有1−2C τ≥0, 即C τ≤
1
, 同时由 2
(
1−C τ−2a 2r 2sin 2σ44
4
)
2
−(1−2C τ)
=4a r sin σ−4(1−C τ) a 2r 2sin 2σ+C 2τ2≤0
中国地质大学(北京)廉海荣编
19
得ar ≤???? (5) (2) 而当C τ>
1
时, 必然有1−C τ−2a 2r 2sin 2σ2
()
2
−(1−2C τ) ≥0, 此时
|λ1,2|=
2.1) 若1−C τ−2a 2r 2sin 2σ≥0, 即a 2r 2≤1−C τ
2
, 此时有
C τ
ρ(G ) = ≤=1
von Neumann条件成立
2.2) 若1−C τ−2a 2r 2sin 2σ≤0, 即a 2r 2≥
1−C τ
2
ρ(G ) =
≤
=−2+4a 2r 2sin 2σ+2C τ 要使von Neumann条件成立, 那么
中国地质大学(北京)廉海荣编
(6) 20
−2+4a 2r 2sin 2σ≤1 (7) 联立(6)-(7), 可得稳定的必要条件???.
5 若矩阵G (σ, r )的两个互异特征值λ1，λ2，那么相应地有两个线性无关的归范特征向
e1 ,e2
e T
1=
iA , 2C τλ1−1]
e T
2=
iA , 2C τλ2−1]
而以e1 ,e2为列向量的行列式之模等于???
则定理4.2.2可知，格式稳定 当ar =1时, 取??? 综上所述，差分格式的稳定性条件是ar
方法二: 类似于抛物型多层差分格式稳定性讨论
将三层格式转换为双层格式组
δk k a 2r 2δk k k +1
t u j =x u j −2τC (u j −u k j )
展开得：u k +1−2u k k −1=a 2r 2(u k k k k +1
j j +u j j +1−2u j +u j −1) −2τC (u j −u k j )
令v k +1
=u k j j , 有：
中国地质大学(北京)廉海荣编 21
⎧⎪u k +1j （
1+2τC ）=a 2r 2u k +22k k
j +1（2+2τC −2a 2r 2）u k j +a r u j −1−v ⎨j ⎪⎩v k +1k
j =u j
令w k j
=⎛⎜u k j
⎞⎜⎝v k ⎟, 写成矩阵形式：
j ⎟⎠
⎛1+2τC 0⎞⎟w k +1⎛2+2τC −2a 2r 2⎜−1⎞j =⎜
⎟w k ⎛a 2r 2
0⎞⎛a 2r 2⎝01⎠⎝10⎠j +⎜⎝00⎟⎠w k
j +1+⎜⎝01 求显格式组的传播因子(二阶矩阵)
w k k e
i σj h
j =q , 可得传播矩阵为 ⎛2+2τC −2a 2r 2(1−cos σh )
G (σ, τ) =⎜1+2−
1⎞
⎜
1+2τC ⎟
⎜τC ⎝
10⎟
⎟⎠
3 求传播矩阵的特征值, 判别von Neumann条件 特征方程：
λ2
−λ
2+2τC −2a 2r 2(1−cos σh ) 1+2τC +1
1+2τC
=0 特征根λ1, 2按模不大于1的充要条件是
2+2τC −2a 2r 2(1−cos σh )
1+2τC
≤2
当ar ≤1 时上式成立，此时von Neumann条件满足
4 利用定理4.2.2给出稳定的充分条件
中国地质大学(北京)廉海荣编 0⎞0⎟⎠
w k
j −1 22
（σ，τ）的两个特征值λ1，λ2当ar
λ1，2=[1+τc −a 2r 2(1−cos σh )]/(1+2τc ) 令A =[1+τc −a 2r 2(1−cos σh )]/(1+2τc ), B =1/(1+2τc )
⎞⎛A −B ⎞⎛00
则(
G-
λ1E )=⎜~⎜⎜1−A −A ⎝1⎝所以对于特征值λ1的特征向量e1=A 同理对于特征值λ2的特征向量e 2e 1 e 2=
A (
=(
A ))
T
T
A 1
=≥δ>0
故由4.4节的稳定性判别定理知差分格式稳定
当ar =1，σ=2N π时
⎛2τC +2
G (σ, τ) =⎜1+2τC
⎜⎜1⎝
−
1⎞
1+2τC ⎟ ⎟⎟0⎠
λ2−λ
2τC +21
+=0
1+2τC 1+2τC
特别地τC →0时时特征方程有特征根：λ1,2=1 此时：G (σ, τ) =⎜
⎛2−1⎞
⎟ 10⎝⎠
⎛k +1−k ⎞
G k (σ, τ) =⎜⎟G (σ, τ) ∞=2k +2→+∞, 无界有定义可知，此时k −k −1⎝⎠
格式不稳定。ar
9. 对一阶常系数方程
中国地质大学(北京)廉海荣编
∂u ∂u +1k +1k
+a =0,试证差分格式αu k j +1+βu j −1=u j ∂t ∂x
23
当α−β≥1时是稳定的;特别地,对格式
1k +1+1k
+1k +1(u j +1+u k j −1) −u j u k j +1−u j −1+a =0 2h τ
当ar ≥1时是稳定的,其中r =τ.
k i σj h
证明：此格式为二层格式. 令u k 代入上式中有， j =v e
v
k +1
e i σjh
v k =i σ(j +1) h
i σ(j −1) h
+βe αe
αe
i σh
则有传播因子：G (σ, τ) =
1
+βe −i σh
由欧拉公式, 即e i σh =cos σh +i sin σh , e −i σh =cos σh −i sin σh , 可得
G (σ, τ) =
1
α(cosσh +i sin σh ) +β(cosσh −i sin σh )
1
=
(α+β) cos σh +i (α−
β)sin σh
显然 G (σ, τ) =, 当α−β≥1时,
G (σ, h ) ≤1，满足冯·诺依曼条件,所以格式是稳定的.
1k +1+1k
+1k +1(u j +1+u k j −1) −u j u k j +1−u j −1k k i σj h +a =0, 有 令u j =v e 代入格式
τ2h
24
中国地质大学(北京)廉海荣编
v
k +1
(m ) =
2e i σjh
v k (m ) [e i σ(j +1) h
+e i σ(j −1) h +ar (e i σ(j +1) h −e i σ(j −1) h
)]
其中网格比r =τ，又由欧拉公式, 可得传播因子为： G (σ, τ) =
1
cos(σh ) +iar sin(σh )
显然： G (σ, τ) =
ar ≥1时,
≥=1 故有G (σ, τ) ≤1, 所以ar ≥1时格式是稳定的.
对于此格式也可以利用第一问的结果来验证, 事实上格式等价于
(1ar +12+2u k 1ar +1+(2−2
) u k +1j −1=u k j j 根据本题的第一部分，可知α=
12+ar 1ar
2, β=2−2
, 此时 −β=
1ar +ar −−ar 2+2−12−ar
2=
2
而当ar ≥1，得α−β=1, 满足α−β≥1，故格式稳定。
中国地质大学(北京)廉海荣编 25