第八章 二維有限差分法數值計算模擬
8.2 邊界條件和有限差分式推導
在此以 TE wave 為例說明，二維模型圖 8-1 上下邊界為電場或是 磁牆，設點方式是距離上下邊界半個格點∆x 2，以此設點方式上下 邊界是電牆或磁牆點會在相同的位置，若是上下邊界電牆，在邊界上 電場內插為零，以上邊界為例u(0,j) =−u(1,j)，若是上下邊界為磁牆，
在邊界上切線電場連續，以上邊界為例u(0,j) =u(1,j)，其中j=1,...,M。
左右方向邊界為吸波邊界條件，解析解形式方程式如式(8-3)到
(8-13)。設點方式如圖 8-4，水平方向設點的作法則是延伸一維不連續
克脈衝函數(Dirac delta function)式(8-21)
)
再將左邊界左邊的點z = −∆z代入式(8-6)，可得一個函數再將不連續
有限差分式，等效介電常數和等效導磁係數以積分的形式求得，如式
圖 8-5 水平方向等效導磁係數
圖 8-6 垂直方向等效導磁效數
和µ 為相對位置向下半個格點、向上半個格點、向右半個格點和向左l 半個格點的導磁係數。I 是元素皆為 1的矩陣大小為 N×M。
d
d ε µ
C 1 1 1
x2
= ∆ (8-47)
u
u ε µ
C 1 1 1
x2
= ∆ (8-48)
r
r ε µ
C 1 1 1 z2
= ∆ (8-49)
l
l ε µ
C 1 1 1 z2
= ∆ (8-50)
µ I µ µ
µ C ε
d u l
r c
2 2 0
2 1 1 )
1 ( 1 )
( 1 1
1 k
x
z +
+
+ ∆
∆ +
= − (8-51)
由於上下邊界是電牆或是磁牆，所以Cu的第一列會反射到Cc的第 1列，電牆u(0,j) =−u(1,j)、磁牆u(0,j) =u(1,j)，其中j=1,...,M，會產生Cu 的第1列為零C (1,1:M)=0u ，同樣的Cd的第N列會反射到Cc第N列，
Cd的第N 列會為零C (N,1:M)=0d 。在左右邊界的部分Cr的第M行會 和右邊界值反射回Cc，同樣的Cl的第1行和左邊界值反射回Cc。
將有限差分式以行設點表示出計算區域，最後的矩陣會成為如此 形式以子矩陣的方式看是腰寛為三的對角矩陣，如式(8-52)。其中的 矩陣是以子矩陣的形式，一行以一個子矩陣表示，例如一行有 N 個 點，會產生子矩陣大小為N×N。
數位置矩陣Cl第1行第j個元素C (j,1)l ，其中j=1,...,N。 對角線元素diag(AL,0)必須還要再加上參數位置矩陣
Cc的第 1行向量C (j,1)c ，還要加上波動方程式中的k02 項。
右對角線元素diag(AL,1)還必須加上參數位置矩陣Cd 第1行元素的N-1個C (1:N-1,1)d 。
左對角線元素diag(AL,-1)還必須加上參數位置矩陣Cu 第1行元素的N-1個C (2:N,1)u 。
A R
若是吸波邊界為密矩陣，若是電牆或是磁牆則是腰寛 為三的對角矩陣。
若是吸波邊界從右邊界條件得到的矩陣T1T⋅P⋅T1第j 列必須要乘上參數位置矩陣Cr第M行第 j個元素
C (j,M)r ，其中j=1,…,N。若是電牆反射參數位置矩陣 -C (1:N,M)r 回C (1:N,M)c ，磁牆反射C (1:N,M)r 回
C (1:N,M)c 。
對角線元素diag(AR,0)必須還要再加上參數位置矩陣 Cc的第 M行向量C (1:N,M)c ，還要加上波動方程式中 的k02項。
右對角線元素diag(AR,1)還必須加上參數位置矩陣Cd 第M行元素的N-1個C (1:N-1,M)d 。
左對角線元素diag(AR,-1)還必須加上參數位置矩陣Cu 第M行元素的N-1個C (2:N,M)u 。
U
D i 為對角線矩陣，i=1,…,M-1。
對角線元素由參數位置矩陣Cr的第 i行表示
C (1:N,i)r 。參數位置矩陣Cr的第M行在右邊界已反射 到AR，所以沒有產生子矩陣DUM。
L
D i
為對角線矩陣，i=2,…,M。
對角線元素由參數位置矩陣Cl的第 i行表示
C (1:N,i)l 。參數位置矩陣Cl的第 1行在左邊界已反射 到AL，所以沒有產生子矩陣D1L。
u i 為行向量大小是N×1，其中i=1,…,M，在計算區域中 設點以行向量為單位，行向量的元素有N個。
Sinc
為入射波產生的項，由左吸波邊界條件產生的一個行 向量大小N×1，若是上下為電場的話，值為純虛數
φ1
⋅
∆
−2jsin(β1i z) ，上下邊界磁場值則為 ψ0
⋅
∆
−2jsin(β1i z) 。
若是以matlab程式的語法表示如下
子矩陣 matlab 程式語法
Ai
i=2,…,M-1
對角線元素diag(Ai,0)=Cc(1:N,i)+k *ones(N,1)02 。 右對角線元素diag(Ai,1)=Cd(1:N-1,i)。
左對角線元素diag(Ai,-1)=C (2:N,i)u 。
AL
A =diag(L C (1:N,1))*l T1T⋅P⋅T1。 對角線元素diag(AL,0)
=diag(AL,0)+Cc(1:N,1)+k *ones(N,1)02 。
右對角線元素diag(AL,1)=diag(AL,1)+C (1:N-1,1)d 。
左對角線元素diag(AL,-1)=diag(AL,-1)+C (2:N,1)u 。
AR
若是吸波邊界
AR=diag(Cr(1:N,M))*T1T⋅P⋅T1 對角線元素diag(AR,0)
=diag(AR,0)+C (1:N,M)+c k *ones(N,1)02 。
右對角線元素diag(AR,1)=diag(AR,1)+Cd(1:N-1,M)。 左對角線元素diag(AR,-1)=diag(AR,-1)+Cu(2:N,M)。 若是電牆或是磁牆。電牆a=-1，磁牆a=1。
對角線元素diag(AR,0)=C (1:N,M)+a*c C (1:N,M)r 。 右對角線元素diag(AR,1)=Cd(1:N-1,M)。
左對角線元素diag(AR,-1)=Cu(2:N,M)。
U
Di i=1,…,M-1
U
Di =C (1:N,i)r 。
L
Di i=2,…,M
L
Di =Cl(1:N,i)。
Sinc 上下邊界電場Sinc=−2jsin(β1i∆z)⋅φ1。 上下邊界磁場Sinc=−2jsin(β1i∆z)⋅ψ0。