第33 节 线性平流方程

1
注釋
第3.3節 線性平流方程
P.164ጣܒҁࢺПแȈөࠉ৯Ϸΰ෬ПਰޟᛧۡܒϷݙ
考慮線性平流方程
0=∂∂+∂∂x
u c t u (1) 其中0>c 。

上式的一個差分方案是
0∆∆1
1=−+−−+x
u u c
t
u u n
j n j n j
n j (2)
(2)式中時間導數項用向前差分方案來近似，空間導數項則用向後差分方案，這個差分方案稱為向前時間上游方案（forward time upsyream scheme ）。

上式可改寫為
n j n j n j u u u )1(11
µµ−+=−+ (3)
其中x t c /∆∆=µ稱為Courant 數。

傳統的數值穩定性分析，大都和流體動力穩定性分析一樣，使用諧波分析法（harmonic analysis ），先假設變數u 具有下面形式的解：
)(ˆ),(t kx i e u
t x u σ−= (4) 其中u
ˆ表示振幅，實數k 表示波數；σ則是複數，它的實部表示頻率，虛部表示成長率。

對離散時空問題來說，
x j x ∆=, t n t ∆=
因而(4)式改寫為
2
)∆∆(ˆt n x j k i n j e u
u σ−= (5) 為方便起見，令t i e ∆σλ−=，此外因為只有一個變數，故可省略振幅u
ˆ，於是(5)式改寫為
x
ikj n n j
e u ∆λ= (6) 對數值穩定性分析來說，使用(6)式所示的解的形式比較方便。

穩定條件可用下式來區分：
1>λ, 數值方案是不穩定的
1=λ, 數值方案是中性的 1<λ, 數值方案是穩定的
不過根據von Karman 穩定性必要條件（見曾1993第81頁），當
)(∆1t O +<λ
時，數值方案也被認為是穩定的。

這個條件允許數值解隨t ∆呈線性增加，但不能更快。

將(6)式所示的解代入(3)式，得到
µµλ−+=−1∆x ik e (7)
上次整理為
x k i x k ∆−∆+−=sin cos 1µµµλ (8)
由上式有
2
2
2222
cos )1(2)1(sin )cos 1(µ
µµµµµµλ+∆−+−=∆+∆+−=x k x
k x k (8)
(8)式改寫為
)cos 1)(1(21cos )1(222122
µµµµµµµλ−−−=∆−++−=x k
3
根據穩定條件1<λ，下式必須成立：
1)cos 1)(1(21<−−−µµµ (9)
由上式有
0)cos 1)(1(2<−−−µµµ (10)
(10)式可再改寫為
0)1(<−−µµ
由此得到穩定條件如下：
10<<µ (11)
這是(3)式的穩定條件。

必須指出，Courant 數x t c /∆∆=µ，故c 必須大於0，在這條件下(3)式所示的方案才是穩定的。

假如0<c ，則(1)式的向前時間上游方案是
0∆∆11=−+−++x
u u c
t
u u n
j
n j n j
n j (12)
上式中時間導數和空間導數項都向前差分來近似。

上式改寫為
n
j n j n j
u x
t c u x t c u )∆∆1(∆∆11
++−=++ (13) (13)式可再改寫為
n j n j n j u u u )1(11
µµ−+=++ (13)
其中x t c /∆∆=µ。

將(6)式所示的解代入(13)式，得到
µµλ−+=1∆x ik e (14)
上次整理為
x k i x k ∆+∆+−=sin cos 1µµµλ (15)
4
用同樣的方法得到
)cos 1)(1(212
µµµλ−−−=
因而穩定條件是
10<<µ (16)
因此，向前時間上游方案，若0>c ，應為(3)式，若0<c ，則應寫為(12)式，這兩個方案的穩定條件為10<<µ，其中x t c /∆∆=µ。

p.164 平流方程的跳步法
考慮線性平流方程
0=+x t u ηη (1)
上式的一個差分方案如下：
0∆2∆21
1
112=−+−−−−+−x u t n j n j n j
n j ηηηη
(2)
上式中時間差分使用跳步法，空間導數則用中差分來近似。

上式可
改寫為
1
,,3,2,,3,2,
0)(1
1112−===−+−−−−+−M j N
n n j n j n j n j L L ηηµηη (3)
這個方案的穩定條件滿足CFL 判據，即Courant 數x t u /∆∆=µ在0到1之間。

差分方案(3)式，跟原來的線性平流方程和向前時間上游方案不一樣，在時間上和空間上都是二階的，故需要兩個初始條件和兩個邊界條件：
M j j j ,,3,210L ===已知,
已知,
ηη (4) N n n
M n ,,1,01L ===已知,
已知,
ηη (5)
(4)式中第1個初始條件是原來平流方程的初始條件，稱為物理初始條件；第2個是按差分方案(37)式進行預報時所需的另一個條件，稱為計算初始條件。

5
若要由觀測資料n j η~決定最佳的初始值，可令下面的目標函數取極小值:
∑∑−=
n j
n j
n j J 2)~(21ηη (6) 必須記得，上式中總和的上下限並不重要，不寫也沒關係。

由於有兩個初始條件，在這裡我們有兩個做法。

第一，直接決定兩個初始
值，也就是說把兩個初始值1j η, 0
j η當做控制變數，此時約束條件是
(3)式。

第二，以(33)式進行預報時所需的計算初始條件1j η以Euler 方案求出，即使用下面的差分方案先求出計算初始條件：
0∆2∆1
1
111
=−+−−−−+−x
u
t
n j n j n j
n j ηηηη, 1=n
即
M j j j j j ,,3,2,
0)(5.00
10101L ==−+−−+ηηµηη (7)
在這情況下，只需物理初始條件就能進行向前預報，也就是說控制變數是0j η，但約束條件除了(33)式外還有(7)式。

首先考慮第一個方法。

將(36)式視為強約束條件, 則需取極小值的目標函數變為
)]([)~(2111112212
2−−−+−=−=−+−−−=∑∑∑∑n j n j n j n j n j N n M j n j n j n j J ηηµηηληη (8) 將上式對m
k η微分，再將上下標(m , k )改回(n , j ) ，得到
1
,,3,2,,1,0),(~11
112−==−−+−−=∂∂+++−+M j N
n J n j n j n j n j n j n j n j
L L λλµλληηη (9)
(9)式將給出伴隨方程和目標函數關於初始值的梯度。

令1,0≠n 時的(9)式等於零，得到下面的伴隨方程：
6
1
,,3,22
,,1,,
~)(11112−=−=−+−−=+++−+M j N N n n
j n j n j n j n j n j L L ηηλλµλλ (10)
當n =N 時，(10)式右邊第1項和第2項出現了2+N j λ, 1
+N j λ，但由(8)
式可知，這兩項並不存在，不過若令
1,,3,2,0,01
2
−===++M j N j
N j
L λλ (11)
則(10)式仍然成立，(11)式是伴隨方程(10)式的終端條件。

另外，當j =2,
j =M 時，(10)式右邊第2項出現了n
1λ, n M 1+λ這兩項，但由(8)式可知，
這一項也不存在，但如果令
N n n M n
,,2,1,0,
0,
01L ===λλ (12)
則(10)式仍然成立，(12)式是伴隨方程的邊界條件。

當n =1和n =0時，(9)式給出目標函數關於兩個初始值的梯度：
1,,3,2),(~212131111−=−−+−−=∂∂+−M j J j j j j j j j L λλµλληηη 1,,3,2),(~111120000−=−−+−−=∂∂+−M j J j j j j j j j L λλµλληηη
由(8)式可知，上式中1j λ和0j λ也是不存在的，故上二式應改寫為
1,,3,2),(~21213111−=−−+−=∂∂+−M j J j j j j j j L λλµληηη (13) 1,,3,2),(~11112000−=−−+−=∂∂+−M j J j j j j j j
L λλµληηη (14)
可以看出，對伴隨方程(10)式進行反向積分時，若再進一步求出1j λ和
0j λ，則這個量正好等於(13)和(14)式的右邊，因此
1,,3,2,,0
01
0−==∂∂=∂∂M j J J j j
j j
L ληλη
(15)
(15a,b)式給出了目標函數關於兩個初始值的梯度。

因此，對伴隨方
程積分時要從n =N 積到n =0，也就是說要計算出1j λ和0j λ，這兩個量正好等於下降算法需用到的目標函數關於兩個初始值的梯度。

7
接著考慮第2個方法。

將(3)式和(7)式當做強約束條件，引進
Lagrange 乘數n j λ和j σ，則目標函數變為
)]
(5.0[)]([)~(210101011
2
1111221
2
2−+−=−−−+−=−=−+−−
−+−−−=∑∑∑∑∑j j j j M j j n j n j n j n j n j N n M j n j n j n j J ηηµηησηηµηηληη (16)
將上式對n j η微分，有
0,111,11112)](5.0[)(~n j j j n j n j n j n j n j n j n j n j
J δσσµσδσλλµλληηη+−+++−+−−+−−−+−−=∂∂ (17) 其中1,n δ和0,n δ是Kronecker 符號。

當0≠n 時，令(47)式等於零，得到
)(~21
21311+−−−+−=j j j j j j λλµληησ )1(=n (18) )(~11
112+++−+−−+−=n j n j n j n j n j n j λλµληηλ , 2,,1,L −=N N n (19) 當0=n 時，由(17)式有 )(5.0~112000+−−−++−=∂∂j j j j j j j J σσµσληηη (20) 上式中已令001==j j λλ，因為它們不存在。

假如伴隨方程(19)式繼續積到1=n ，即
)(~11
112+++−+−−+−=n j n j n j n j n j n j λλµληηλ , 1,2,,1,L −=N N n (21) 則我們發現1j j λσ=，於是(20)式可改寫為 )(5.0~111112000+−−−++−=∂∂j j j j j j j
J λλµλληηη (22) (22)式雖可再簡化，但不能寫成像(15)式那樣簡單的形式。