擬譜法 (Pseudo-spectral methods)[ 1] 也稱為離散變數表示 (discrete variable representation、DVR)法,是在應用數學 及計算科學 中求解偏微分方程 用的數值分析 方法。擬譜法和譜方法 有密切關係,但在譜方法中基底 函數中使用了擬譜的基底函數,也就是可以在分割網格上表示的函數。此作法簡化一些運算子的計算,在使用快速演算法(例如快速傅立葉變換 )時可以加速計算速度。
說明
考慮以下的初值問題
i
∂
∂
t
ψ
(
x
,
t
)
=
[
−
∂
2
∂
x
2
+
V
(
x
)
]
ψ
(
x
,
t
)
,
ψ
(
t
0
)
=
ψ
0
{\displaystyle i{\frac {\partial }{\partial t}}\psi (x,t)={\Bigl [}-{\frac {\partial ^{2}}{\partial x^{2}}}+V(x){\Bigr ]}\psi (x,t),\qquad \qquad \psi (t_{0})=\psi _{0}}
有週期性條件
ψ
(
x
+
1
,
t
)
=
ψ
(
x
,
t
)
{\displaystyle \psi (x+1,t)=\psi (x,t)}
。這個例子是勢能為
V
(
x
)
{\displaystyle V(x)}
之粒子的薛丁格方程 ,不過其結構可適用到其他應用。在許多實務上的偏微分方程中,其中有一項是和導數(例如是動能相關的量)有關,另一項則是和另一個函數(此處為勢能)的乘積。
在譜方法中,其解
ψ
{\displaystyle \psi }
會展開為一組適合基底函數的組合,例如平面波。
ψ
(
x
,
t
)
=
1
2
π
∑
n
c
n
(
t
)
e
2
π
i
n
x
.
{\displaystyle \psi (x,t)={\frac {1}{\sqrt {2\pi }}}\sum _{n}c_{n}(t)e^{2\pi inx}.}
將解代入,並且計算係數的方程,可以得到係數的常微分方程
i
d
d
t
c
n
(
t
)
=
(
2
π
n
)
2
c
n
+
∑
k
V
n
−
k
c
k
,
{\displaystyle i{\frac {d}{dt}}c_{n}(t)=(2\pi n)^{2}c_{n}+\sum _{k}V_{n-k}c_{k},}
而其元素
V
n
k
{\displaystyle V_{nk}}
可以透過顯式的傅立葉轉換求得
V
n
−
k
=
∫
0
1
V
(
x
)
e
2
π
i
(
k
−
n
)
x
d
x
.
{\displaystyle V_{n-k}=\int _{0}^{1}V(x)\ e^{2\pi i(k-n)x}dx.}
若將
N
{\displaystyle N}
基底函數的展開到一定項次後截斷,並且找
c
n
(
t
)
{\displaystyle c_{n}(t)}
的解,就可以得到偏微分方程的解。一般而言會用數值方法 進行,例如龍格-庫塔法 。為了數值解,常微分方程的右側需重覆計算在許多不同時間間隔下的值。此時,譜方法有個和勢能項
V
(
x
)
{\displaystyle V(x)}
有關的主要問題。
在譜方法下,和勢能函數
V
(
x
)
{\displaystyle V(x)}
的相乘會轉換為向量和矩陣的乘法,其複雜度是
N
2
{\displaystyle N^{2}}
,而且在求解係數的微分方程時,需要另外去計算矩陣元素
V
n
k
{\displaystyle V_{nk}}
,這也需要時間。
在擬譜法中,會用不同的方式來計算。給定係數
c
n
(
t
)
{\displaystyle c_{n}(t)}
,會用反離散傅立葉轉換來計算函數
ψ
{\displaystyle \psi }
在離散格點
x
j
=
2
π
j
/
N
{\displaystyle x_{j}=2\pi j/N}
下的值。在格點上,計算函數的乘積
ψ
′
(
x
i
,
t
)
=
V
(
x
i
)
ψ
(
x
i
,
t
)
{\displaystyle \psi '(x_{i},t)=V(x_{i})\psi (x_{i},t)}
,再用傅立葉轉換轉換回來,可以得到一組新的係數
c
n
′
(
t
)
{\displaystyle c'_{n}(t)}
,來代替矩陣乘積運算
∑
k
V
n
−
k
c
k
(
t
)
{\displaystyle \sum _{k}V_{n-k}c_{k}(t)}
。
可以證明二個方法有類似的精準度,而且擬譜法可以使用快速傅立葉轉換,其時間複雜度為
O
(
N
ln
N
)
{\displaystyle O(N\ln N)}
,理論上比矩陣乘法要快很多。而且,可以直接計算函數
V
(
x
)
{\displaystyle V(x)}
,不用再經過額外的積分運算。
技術討論
若用比較抽象的方式來描述,擬譜法是處理在偏微分方程中二個函數
V
(
x
)
{\displaystyle V(x)}
和
f
(
x
)
{\displaystyle f(x)}
的乘積。為了簡化表示式,省略掉函數中時間的自變數。在概念上,擬譜法包括三個步驟:
將
f
(
x
)
{\displaystyle f(x)}
以及
f
~
(
x
)
=
V
(
x
)
f
(
x
)
{\displaystyle {\tilde {f}}(x)=V(x)f(x)}
擴展為由有限多個基底函數形成的組合(此為譜方法 )。
針對一組給定的基底函數,找到一組分割方式可以將基底函數的純量積轉換為在格點上的加權和。
在每一個格點將
V
{\displaystyle V}
和
f
{\displaystyle f}
相乘,來計算函數的乘積。
基底的展開
函數
f
{\displaystyle f}
和
f
~
{\displaystyle {\tilde {f}}}
可以用一組有限基底的函數來擴展
{
ϕ
n
}
n
=
0
,
…
,
N
{\displaystyle \{\phi _{n}\}_{n=0,\ldots ,N}}
:
f
(
x
)
=
∑
n
=
0
N
c
n
ϕ
n
(
x
)
{\displaystyle f(x)=\sum _{n=0}^{N}c_{n}\phi _{n}(x)}
f
~
(
x
)
=
∑
n
=
0
N
c
~
n
ϕ
n
(
x
)
{\displaystyle {\tilde {f}}(x)=\sum _{n=0}^{N}{\tilde {c}}_{n}\phi _{n}(x)}
為了簡化起見,令基底是正交且正規化的,
⟨
ϕ
n
,
ϕ
m
⟩
=
δ
n
m
{\displaystyle \langle \phi _{n},\phi _{m}\rangle =\delta _{nm}}
,利用內積
⟨
f
,
g
⟩
=
∫
a
b
f
(
x
)
g
(
x
)
¯
d
x
{\displaystyle \langle f,g\rangle =\int _{a}^{b}f(x){\overline {g(x)}}dx}
配合適當的邊界
a
,
b
{\displaystyle a,b}
,其係數為
c
n
=
⟨
f
,
ϕ
n
⟩
{\displaystyle c_{n}=\langle f,\phi _{n}\rangle }
c
~
n
=
⟨
f
~
,
ϕ
n
⟩
{\displaystyle {\tilde {c}}_{n}=\langle {\tilde {f}},\phi _{n}\rangle }
配合一些計算可得
c
~
n
=
∑
m
=
0
N
V
n
−
m
c
m
{\displaystyle {\tilde {c}}_{n}=\sum _{m=0}^{N}V_{n-m}c_{m}}
而
V
n
−
m
=
⟨
V
ϕ
m
,
ϕ
n
⟩
{\displaystyle V_{n-m}=\langle V\phi _{m},\phi _{n}\rangle }
。這是譜方法的基礎。為了區分
ϕ
n
{\displaystyle \phi _{n}}
的基底以及正交的基底,有時會將上述擴展稱為有限基底表示(Finite Basis Representation、FBR)。
分割
針對給定基底
{
ϕ
n
}
{\displaystyle \{\phi _{n}\}}
以及
N
+
1
{\displaystyle N+1}
個基底函數,可以設法找到分割方式,也就是
N
+
1
{\displaystyle N+1}
個點以及加權,使得
⟨
ϕ
n
,
ϕ
m
⟩
=
∑
i
=
0
N
w
i
ϕ
n
(
x
i
)
ϕ
m
(
x
i
)
¯
n
,
m
=
0
,
…
,
N
{\displaystyle \langle \phi _{n},\phi _{m}\rangle =\sum _{i=0}^{N}w_{i}\phi _{n}(x_{i}){\overline {\phi _{m}(x_{i})}}\qquad \qquad n,m=0,\ldots ,N}
特別的例子包括多項式的高斯求積 以及平面波的離散傅立葉變換 ,特別需注意的是格點及加權
x
i
,
w
i
{\displaystyle x_{i},w_{i}}
都是基底及數量
N
{\displaystyle N}
的函數。
利用分割方式,可以透過格點上的值,以另一種方式來表示函數
f
(
x
)
,
f
~
(
x
)
{\displaystyle f(x),{\tilde {f}}(x)}
的數值。此表示法有時稱為離散變數表示法(Discrete Variable Representation、DVR),完全等效於基底的展開。
f
(
x
i
)
=
∑
n
=
0
N
c
n
ϕ
n
(
x
i
)
{\displaystyle f(x_{i})=\sum _{n=0}^{N}c_{n}\phi _{n}(x_{i})}
c
n
=
⟨
f
,
ϕ
n
⟩
=
∑
n
=
0
N
w
i
f
(
x
i
)
ϕ
n
(
x
i
)
¯
{\displaystyle c_{n}=\langle f,\phi _{n}\rangle =\sum _{n=0}^{N}w_{i}f(x_{i}){\overline {\phi _{n}(x_{i})}}}
相乘
和函數
V
(
x
)
{\displaystyle V(x)}
的相乘會在格點上進行
f
~
(
x
i
)
=
V
(
x
i
)
f
(
x
i
)
.
{\displaystyle {\tilde {f}}(x_{i})=V(x_{i})f(x_{i}).}
一般來說這裡會有一些近似,可以計算其中一個係數
c
~
n
{\displaystyle {\tilde {c}}_{n}}
:
c
~
n
=
⟨
f
~
,
ϕ
n
⟩
=
∑
i
w
i
f
~
(
x
i
)
ϕ
n
(
x
i
)
¯
=
∑
i
w
i
V
(
x
i
)
f
(
x
i
)
ϕ
n
(
x
i
)
¯
{\displaystyle {\tilde {c}}_{n}=\langle {\tilde {f}},\phi _{n}\rangle =\sum _{i}w_{i}{\tilde {f}}(x_{i}){\overline {\phi _{n}(x_{i})}}=\sum _{i}w_{i}V(x_{i})f(x_{i}){\overline {\phi _{n}(x_{i})}}}
利用譜方法,對應的係數會是
c
~
n
=
⟨
V
f
,
ϕ
n
⟩
{\displaystyle {\tilde {c}}_{n}=\langle Vf,\phi _{n}\rangle }
。擬譜法則會用以下皂的近似來處理
⟨
V
f
,
ϕ
n
⟩
≈
∑
i
w
i
V
(
x
i
)
f
(
x
i
)
ϕ
n
(
x
i
)
¯
.
{\displaystyle \langle Vf,\phi _{n}\rangle \approx \sum _{i}w_{i}V(x_{i})f(x_{i}){\overline {\phi _{n}(x_{i})}}.}
若乘積可以用
V
f
{\displaystyle Vf}
給定的有限基底函數組合來表現,則上式在給定分割方式上會完全正確。
特殊的擬譜架構
傅立葉法
若問題中有週期性的邊界條件,其週期為
[
0
,
L
]
{\displaystyle [0,L]}
,基底函數可以用平面波來產生,
ϕ
n
(
x
)
=
1
L
e
−
ı
k
n
x
{\displaystyle \phi _{n}(x)={\frac {1}{\sqrt {L}}}e^{-\imath k_{n}x}}
其中
k
n
=
(
−
1
)
n
⌈
n
/
2
⌉
2
π
/
L
{\displaystyle k_{n}=(-1)^{n}\lceil n/2\rceil 2\pi /L}
,而
⌈
⋅
⌉
{\displaystyle \lceil \cdot \rceil }
是取整函數 。
在
n
max
=
N
{\displaystyle n_{\text{max}}=N}
處截斷的分割為離散傅立葉轉換 ,格點會平均分佈
x
i
=
i
Δ
x
{\displaystyle x_{i}=i\Delta x}
,間隔為
Δ
x
=
L
/
(
N
+
1
)
{\displaystyle \Delta x=L/(N+1)}
,各點的加權會是相同旳定值
w
i
=
Δ
x
{\displaystyle w_{i}=\Delta x}
。
在討論誤差時,需注意到平面波的乘積也是平面波,
ϕ
a
+
ϕ
b
=
ϕ
c
{\displaystyle \phi _{a}+\phi _{b}=\phi _{c}}
,
c
≤
a
+
b
{\displaystyle c\leq a+b}
。因此,定量來說,若函數
f
(
x
)
,
V
(
x
)
{\displaystyle f(x),V(x)}
可以用
N
f
,
N
V
{\displaystyle N_{f},N_{V}}
基底函數足夠準確的呈現,只要用
N
f
+
N
V
{\displaystyle N_{f}+N_{V}}
個基底函數,即可用擬譜法得到足夠準確的結果。
平面波的擴展一般效果較差,需要許多的基底函數才能收斂。不過。基底展開和格點表示的轉換可以用快速傅立葉轉換 進行,其時間複雜度較低,為
N
ln
N
{\displaystyle N\ln N}
。因此,平面波是擬譜法中常用的一種基底函數。
多項式
另一種常見的展開方式是多項式,此處會使用高斯求積 (Gaussian quadrature),其中提到可以找到加權係數
w
i
{\displaystyle w_{i}}
及格點
x
i
{\displaystyle x_{i}}
使得
∫
a
b
w
(
x
)
p
(
x
)
d
x
=
∑
i
=
0
N
w
i
p
(
x
i
)
{\displaystyle \int _{a}^{b}w(x)p(x)dx=\sum _{i=0}^{N}w_{i}p(x_{i})}
對任意的
2
N
+
1
{\displaystyle 2N+1}
次或是更低次的多項式
p
(
x
)
{\displaystyle p(x)}
都成立。一般而言,加權函數
w
(
x
)
{\displaystyle w(x)}
及範圍
a
,
b
{\displaystyle a,b}
都是根據特定問題所選定的,因此會選擇幾種分割方式中的一種。若要用在擬譜法,需選擇基底函數為
ϕ
n
(
x
)
=
w
(
x
)
P
n
(
x
)
{\displaystyle \phi _{n}(x)={\sqrt {w(x)}}P_{n}(x)}
,其中
P
n
{\displaystyle P_{n}}
為
n
{\displaystyle n}
階多項式,有以下特性
∫
a
b
w
(
x
)
P
n
(
x
)
P
m
(
x
)
d
x
=
δ
m
n
.
{\displaystyle \int _{a}^{b}w(x)P_{n}(x)P_{m}(x)dx=\delta _{mn}.}
在上述條件下,the
ϕ
n
{\displaystyle \phi _{n}}
會是正交基底,其內積為
⟨
f
,
g
⟩
=
∫
a
b
f
(
x
)
g
(
x
)
¯
d
x
{\displaystyle \langle f,g\rangle =\int _{a}^{b}f(x){\overline {g(x)}}dx}
。此基底以及分割點可以用在擬譜法中。
有關其誤差,若
f
{\displaystyle f}
可以用
N
f
{\displaystyle N_{f}}
個基底函數很好的呈現,而
V
{\displaystyle V}
可以用
N
V
{\displaystyle N_{V}}
階多項式很好的呈現,則其積可以用前
N
f
+
N
V
{\displaystyle N_{f}+N_{V}}
個基底函數很好的呈現。且擬譜法用該數量的基底函數,會有足夠準確的結果。
在一些標準問題中,會出現這些多項式。例如量子簡諧振動子可以擴展為埃爾米特多項式 ,而在轉動問題中,會用雅可比多項式 來定義相關的勒壤得多項式 。
參考資料
Orszag, Steven A. Numerical Methods for the Simulation of Turbulence. Physics of Fluids. 1969, 12 (12): II–250. doi:10.1063/1.1692445 .
Gottlieb, David; Orszag, Steven A. Numerical analysis of spectral methods : theory and applications 5. print. Philadelphia, Pa.: Society for Industrial and Applied Mathematics. 1989. ISBN 978-0898710236 .
Hesthaven, Jan S.; Gottlieb, Sigal; Gottlieb, David. Spectral methods for time-dependent problems 1. publ. Cambridge [u.a.]: Cambridge Univ. Press. 2007. ISBN 9780521792110 .
Trefethen, Lloyd N. Spectral methods in MATLAB 3rd. repr. Philadelphia, Pa: SIAM. 2000. ISBN 978-0-89871-465-4 .
Fornberg, Bengt. A Practical Guide to Pseudospectral Methods. Cambridge: Cambridge University Press. 1996. ISBN 9780511626357 .
Boyd, John P. Chebyshev and Fourier spectral methods 2nd ed., rev. Mineola, N.Y.: Dover Publications. 2001 [2019-07-17 ] . ISBN 978-0486411835 . (原始內容存檔 於2019-06-24).
Funaro, Daniele. Polynomial approximation of differential equations . Berlin: Springer-Verlag. 1992 [2019-07-17 ] . ISBN 978-3-540-46783-0 . (原始內容存檔 於2011-07-22).
de Frutos, Javier; Novo, Julia. A Spectral Element Method for the Navier--Stokes Equations with Improved Accuracy. SIAM Journal on Numerical Analysis. January 2000, 38 (3): 799–819. doi:10.1137/S0036142999351984 .
Claudio, Canuto; M. Yousuff, Hussaini; Alfio, Quarteroni; Thomas A., Zang. Spectral methods fundamentals in single domains. Berlin: Springer-Verlag. 2006. ISBN 978-3-540-30726-6 .
Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP. Section 20.7. Spectral Methods . Numerical Recipes: The Art of Scientific Computing 3rd. New York: Cambridge University Press. 2007 [2019-07-17 ] . ISBN 978-0-521-88068-8 . (原始內容存檔 於2011-08-11).