偽譜法基本原理

偽譜法適用條件：如果你想在一個簡單的領域上高精度地解決一個ODE或PDE，且如果定義問題的數據是平滑的，那么偽譜法通常是最好的工具，偽譜法通常可以達到10位數的精度，其中有限差分或有限元方法將得到2到3位數。在較低的精度下，偽譜法比其他方法需要更少的計算機內存。

在偽譜法中，偏微分方程的解是由一組正交函數$\{{\phi }_j(x){\}}_{j=0}^{\infty }$的有限線性組合來近似的:
$$\begin{array}{c}{f}_N(x)=\sum _{k=0}^N{\hat{f}}_k{\phi }_k(x)\end{array}\text{\hspace{1em}?}$$其中${\hat{f}}_k$被稱為$f$相對于基函數$\{{\phi }_k\}$的譜系數，${\hat{f}}_k$是$f$在${\phi }_k$“方向”上的度量，可以通過一個積分來計算
$$\begin{array}{c}\begin{array}{r}{\hat{f}}_k=\int f(x){\phi }_k(x)\omega (x)\mathrm{d}x\end{array}\end{array}\text{\hspace{1em}?}$$其中，$\omega (x)$是一個與基函數$\{{\phi }_k\}$相關聯的權重函數。
假設，我們考慮的是在$[0,2\pi ]$內的周期函數，其中一個基函數是由正弦和余弦組成的傅里葉基函數，如$\{\cos kx，\sin kx{\}}_{k=0}^{\infty }$，其簡潔的復數形式為$\{e^{ikx}{\}}_{k=?\infty }^{\infty }$，如果$f$滿足一些條件，如連續性和有界可微，可以證明由$f_N(x)={\sum }_{k=0}^N{\hat{f}}_k{e}^{ikx}$其中${\hat{f}}_k=\int f(x)e^{ikx}\mathrm{d}x$給出的,$f$隨著$N$的增加一致收斂于$f$。
$f$的導數近似為$f_N^{{}^{\prime }}$，由下式給出：
$$\begin{array}{c}\begin{array}{r}{f}_N^{{}^{\prime }}(x)=\sum _{k=0}^N{\hat{f}}_k{\phi }_k^{{}^{\prime }}(x)\end{array}\end{array}\text{\hspace{1em}?}$$$f_N^{{}^{\prime }}$被稱為$f$的伽遼金譜導數。一旦我們知道了導數${\phi }_k^{{}^{\prime }}$（一般的遞推公式）的細節，我們就可以計算出公式$(3)$，并找出系數${\hat{f}}_k$來導出$f_N^{{}^{\prime }}(x)={\sum }_{k=0}^N{\hat{f}}_k^{{}^{\prime }}{\phi }_k(x)$。
例如，設$f$是$[0,2\pi ]$中的一個偶數周期函數。近似$f$的一個合適的基函數是傅里葉基函數。如果我們寫的話：
$$\begin{array}{c}\begin{array}{r}{f}_N(x)=\sum _{k=0}^N{\hat{f}}_k{e}^{ikx}\end{array}\end{array}\text{\hspace{1em}?}$$我們很容易看到，$f$的一階和二階導數的近似是：
$$\begin{array}{c}\begin{array}{r}{f}_N^{{}^{\prime }}(x)=\sum _{k=0}^N\left(ik{\hat{f}}_k\right)e^{ikx}\\ f_N^{{}^{\prime \prime }}(x)=\sum _{k=0}^N\left(?k^2{\hat{f}}_k\right)e^{ikx}\end{array}\end{array}\text{\hspace{1em}?}$$這就是說，一旦我們知道了譜系數${\hat{f}}_k$，我們就通過將每個${\hat{f}}_k$乘以$?k^2$得到了$f^{{}^{\prime \prime }}$的近似值。例如使用偽譜法，將譜系數分別乘以$ik,(ik{)}^2,(ik{)}^3$來求解函數$f(x)=\sin (\pi (x+1))e^{\cos (\pi (x+1))}$的一次、二次與三次導數的結果如下：

偽譜法求導結果

偽譜法求解聲波方程

聲波方程為：
$$\begin{array}{c}\begin{array}{r}{?}_t^{2}u=v^2{?}^{2}u=v^2({?}_x^2+{?}_y^2)u\end{array}\end{array}\text{\hspace{1em}?}$$將上式對空間的偏導數用偽譜法表示可得：
$$\begin{array}{c}\begin{array}{r}{v}^2({?}_x^2+{?}_y^2)u=v^2{F}^{?1}\left[?(k_x^2+k_y^2)F(u)\right]\end{array}\end{array}\text{\hspace{1em}?}$$其中$F$為傅氏正變換，$F^{?1}$為傅氏逆變換。對時間的偏導數項用二階中心差分表示為：
$$\begin{array}{c}\begin{array}{r}{?}_t^{2}u=\frac{u^{n+1}?2u^n+u^{n?1}}{\Delta t^2}\end{array}\end{array}\text{\hspace{1em}?}$$將（1）和（8）式代入到（6）中去，可得：
$$\begin{array}{c}\begin{array}{rl}{?}_t^{2}u& =\frac{u^{n+1}?2u^n+u^{n?1}}{\Delta t^2}=v^2{F}^{?1}\left[?(k_x^2+k_y^2)F(u)\right]\\ u^n& =(\Delta t^2{v}^2)F^{?1}\left[?(k_x^2+k_y^2)F(u)\right]+2u^n?u^{n?1}\end{array}\end{array}\text{\hspace{1em}?}$$

傅里葉變換討論

傅里葉變換性質

1.物理空間的離散對應傅里葉空間的有界。
2.物理空間的有界對應傅里葉空間的離散。
對于結論一，我們假定兩個復函數為$f(x)=e^{i{k}_{1}x},g(x)=e^{i{k}_{2}x}$，且其中的$k_1\ne k_2$，$x$是定義在離散空間$h\mathbf{Z}$的網格點$x_j=jh,j\in \mathbf{Z}$上的，則：
$$\begin{array}{c}\begin{array}{r}{f}_j=e^{i{k}_1{x}_j}\\ g_j=e^{i{k}_2{x}_j}\end{array}\end{array}\text{\hspace{1em}?}$$如果$k_1?k_2$是$\frac{2\pi }{h}$的倍數，則對于每個$j$有$f_j=g_j$。由此可見，對于任何復指數$e^{ikx}$，有無限的其他復指數與它匹配的網格$h\mathbf{Z}$ 。因此，在長度為$\frac{2\pi }{h}$的區間內測量網格的波數就足夠了，由于對稱性的原因，我們選擇區間$[?\frac{\pi }{h},\frac{\pi }{h}]$。
對于結論二，可以通過傅里葉正反變換的對稱性來理解，因為在傅里葉變換中，時間和空間是等價的。在傅里葉空間的離散，代入到傅氏反變換中，也是需要計算復指數函數的，同理也會導致物理空間的有界$[0,2\pi ]$。
假定空間上的網格點數$N$是偶數，則網格點的間距是$h=2\pi /N$，可得：
$$\begin{array}{c}\begin{array}{r}\frac{\pi }{h}=\frac{N}{2}\end{array}\end{array}\text{\hspace{1em}?}$$現在，傅里葉域是離散的和有界的。這是因為物理空間中的波必須在區間$[0,2\pi ]$上是周期性的，并且只有具有整數波數的$e^{ikx}$具有所需的周期$2\pi$。

aliasing phenomenon

快速傅里葉變換（fft）的性質

這里討論的是MATLAB中使用的fft和ifft。在使用這些程序時必須稍微小心一點，因為MATLAB假設波數的順序與我們的不同。在我們的符號中使用的波數表示為：
$$\begin{array}{c}\begin{array}{r}{\hat{f}}_{?\frac{N}{2}+1},{\hat{f}}_{?\frac{N}{2}+2}?,{\hat{f}}_{\frac{N}{2}}\end{array}\end{array}\text{\hspace{1em}?}$$因為fft的使用使用了一個復雜的轉換用來匹配該變換滿足對稱特性${\hat{f}}_j={\hat{f}}_{?j}$,這是由大多數傅里葉變換所利用的，它只存儲了大約一半的數組，減少了保持系數所需的內存，并使速度提高了大約兩倍，此時傅里葉空間中的傅里葉系數為：
$$\begin{array}{c}\begin{array}{r}{\hat{f}}_0,{\hat{f}}_1,?,{\hat{f}}_{\frac{N}{2}}\end{array}\end{array}\text{\hspace{1em}?}$$

傅里葉系數排列

設$L=Nh$表示仿真域的物理長度。然后，真正的波數數組就采用了這種形式
$$\begin{array}{c}\begin{array}{r}k(i)=?\begin{array}{ll}\frac{2\pi }{L}i,& i=1,?,\frac{N}{2}\\ \frac{2\pi }{L}(?N+i)& i=\frac{N}{2}+1,?,N\end{array}\end{array}\end{array}\text{\hspace{1em}?}$$下圖是16點，間距為1的傅里葉變換系數的分布情況：

$${\hat{f}}_k$$

偽譜法的穩定性

當用偽譜法數值求解時變偏微分方程時，其模式通常是相同的：空間上的譜微分，時間上的有限差分。例如，人們可以通過歐拉格式來執行時間步進。
以變系數波動方程為例：
$$\begin{array}{c}\begin{array}{r}{?}_{t}u=c{?}_{x}u,c={\sin }^2(x?1)+0.2,u_0=e^{?100(x?1{)}^2}\end{array}\end{array}\text{\hspace{1em}?}$$求解它的迭代格式為：
$$\begin{array}{c}\begin{array}{r}{u}^{n+1}=u^n?2\Delta tc{F}^{?1}[ikF(u)]\end{array}\end{array}\text{\hspace{1em}?}$$當時間步長為$\frac{h}{4}=0.0491$時，它是穩定的。

穩定情況

但是當$\Delta t>\frac{1.9}{N}$時，它的解就變得不穩定了。這里的不穩定性是由于時間離散引起的，可以通過適當的隱式離散化來規避。

不穩定情況

偽譜法的頻散現象

與有限差分法相比，傅里葉方法最顯著的特點之一是數值色散小。設置同樣的參數進行偽譜法與三點、五點有限差分實驗，通過對比波傳播相同時間到達相同位置情況下頻率不同引起的頻散現象來研究各種算法的頻散特性，下面顯示快照允許我們得到兩個重要結論，實驗結果如下：

從圖中可以看出隨著子波主頻的不斷增加，不論是三點還是五點有限差分，它們的頻散現象愈發嚴重，且五點有限差分頻散較三點有所改善。相比之下，偽譜法的頻散現象隨主頻增大不是很嚴重，幾乎不存在頻散。偽譜法理論上不會出現頻散現象，我們在數值實驗中出現頻散是因為時間上作了差分離散導致的頻散，跟穩定性一樣可以使用其他形式的時間離散來進行改進。
從c組圖中可以看出，有限差分解具有明顯的強各向異性色散行為。最準確的方向出現在$\theta =45°$，偽譜法解沒有表現出明顯的色散，但最重要的是，它似乎沒有方向依賴性。也就是說，誤差是各向同性的。
關于頻散性的結論，下圖也可以明顯說明：

頻散性比較

吉布斯現象

吉布斯現象是由亨利·威爾伯拉罕在1848年發現的，然后由j·威拉德·吉布斯在1899年重新發現。
對于具有不連續的周期信號，如果通過加傅立葉級數來重構信號，則在邊緣附近會出現超調。這些超調以遠離邊緣的阻尼振蕩方式向外衰減。這被稱為GIBBS現象，如下圖所示。

因為，當$0<x<\pi$時$f(x)=1$，當$\pi <x<2\pi$時$f(x)=?1$(在$(0,2\pi )$范圍外重復)定義的方波具有傅立葉級數展開:
$$\begin{array}{c}f(x)=\frac{4}{\pi }\sum _{n=1,3,5,?}\frac{1}{n}\sin (nx)\end{array}\text{\hspace{1em}?}$$更一般的,隨著N增加，部分起伏就向不連續點壓縮，但是對任何有限的N值，起伏的峰值大小保持不變，以函數$f(x)=x$為例，其周期為$2\pi$，某個定義區間為$[?\pi ,\pi ]$,其圖像為：

且它的傅里葉展開級數的系數為：
$$\begin{array}{c}\begin{array}{rl}{b}_n& =(\frac{1}{\pi }){\int }_{?\pi }^{\pi }x\sin (nx)\mathrm{d}x\\ & =(?1{)}^{n+1}(\frac{2}{n})\end{array}\end{array}\text{\hspace{1em}?}$$則它的傅里葉級數近似結果如下：

其誤差表示為：

從圖中可以看出其誤差最大值保持不變，約為跳變值的0.09倍。
當使用偽譜法求導數是使用的不是函數一整個周期時，求導結果也會出現吉布斯現象，例如在區間$[0,2\pi ]$上求取$\cos (x)$的導數（它的一個完整周期應該是$[0,2\pi )$，不包含$2\pi$）時會出現吉布斯現象。

偽譜法對sin函數求導

總結

以上是生活随笔為你收集整理的伪谱法求解PDE的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。