一、$\tt Toeplitz$矩陣與循環（$\tt Circulant$）矩陣
定義
為$n\times n$階循環矩陣。
定義 $T_n(i,j)=t_{j-i} $ 為$n\times n$ 階$\tt Toeplitz$矩陣
通過令矩陣$B_n＝$
從而構造出$2n\times 2n$階循環矩陣
假設有一$n\times 1$階列向量$\bf u$
其中，$C_{2n}$可以由快速傅里葉對角化
其中$\bf c$表示$C_{2n}$矩陣的第一列元素，$\bf F$ 表示快速傅里葉($\tt fft$)變換，$\bf F^{-1}$ 表示快速傅里葉($\tt ifft$)逆變換。進一步可寫成
因此，計算$\bf T_n u$等價於計算
查閱文獻我們知道，直接計算$\bf T_n u$的存儲量和計算量分別為$O(n^2)$和$O(n^3)$,但是利用快速傅里葉計算可以將存儲量和計算量分別降為$O(n)$和$O(n \log_2 n)$.
從以下數據可以更直觀的看出FFT顯著的優點
二、數值應用
- 考慮一維橢圓方程
$$－\Delta u=f(x),\qquad a<x<b,\tag{1}$$
滿足齊次$Dirichlet$邊界條件。
對$x\in [a,b]$一致網格剖分：$a=x_0<x_1,\cdots,x_M=b$,$h=\frac{b-a}{M}$。$U,u$分別表示數值解和真解。應用二階中心差分逼近二階導數
$$\Delta u(x_i)= \frac{u(x_{i-1})-2u(x_i)+u(x_{i+1}) }{h^2}+O(h^2).\tag{2}$$
由(2)式可得求解方程(1)的數值格式的矩陣形式
$$A{\bf U}＝\widehat{f}.\tag{3}$$
其中
$$A=\tt －\frac{1}{h^2}toeplitz([-2,1,zeros(1,M-3)]),$$
$$\widehat{f}=( f(x_1),f(x_2),\cdots,f(x_{M-1}) )^T.$$
$${\bf U}=(u_1,u_2,\cdots,u_{M-1})^T.$$
- 考慮二維橢圓方程
$$－\Delta u=f({\bf x,y}),\qquad {(\bf x,y)}\in (a,b)\times (c,d),\tag{4}$$
對$x\in [a,b]$一致網格剖分：$a=x_0<x_1,\cdots,<x_{M_1}=b$,$h_1=\frac{b-a}{M_1}$，$c=y_0<y_1,\cdots,<y_{M_2}=d$,$h_2=\frac{d-c}{M_2}$。$U,u$分別表示數值解和真解。應用二階中心差分逼近二階導數
$$\Delta u(x_i,y_j)= \frac{u(x_{i-1},y_j)-2u(x_i,y_j)+u(x_{i+1},y_j) }{h_1^2}+ \frac{u(x_i,y_{j-1})-2u(x_i,y_j)+u(x_i,y_{j+1}) }{h_2^2}+O(h_1^2+h_2^2).\tag{5}$$
由(5)式可得求解方程(4)的數值格式的矩陣形式
$$A{\bf U}＝\widehat{f}.\tag{6}$$
其中
$$A_1=\tt toeplitz([-2,1,zeros(1,M_1-3)]),$$
$$A_2=\tt toeplitz([-2,1,zeros(1,M_2-3)]),$$
$$ A_x = -\tt \frac{1}{h_1^2} I_{M_2-1} \bigotimes A_1 ,\mbox{(非toeplitz矩陣)}$$
注意到：
$$ I_{M_2-1} \bigotimes A_1U = reshape\Big( A_1 reshape( U,M_1-1,M_2-1 ),( M_1-1 )(M_2-1),1 \Big). $$
$$ A_y = -\tt \frac{1}{h_2^2} A_2 \bigotimes I_{M_1-1} ,$$
$$A = A_x+A_y,$$
$$\widehat{f}=( f(x_1,y_1),f(x_2,y_1),\cdots,f(x_{M_1-1},y_1) , f(x_1,y_2),f(x_2,y_2),\cdots,f(x_{M_1-1},y_2), \cdots\cdots, f(x_1,y_{M_2-1}),f(x_2,,y_{M_2-1}),\cdots,f(x_{M_1-1},,y_{M_2-1}) )^T.$$
$${\bf U}=(u_{1,1},u_{2,1},\cdots,u_{M_1-1,1},u_{1,2},u_{2,2},\cdots,u_{M_1-1,2},\cdots\cdots,u_{1,M_2-1},u_{2,M_2-1},\cdots,u_{M_1-1,M_2-1})^T.$$
由數值格式（3），（6）顯然可知，當空間網格剖分很大時，矩陣乘法的計算量將會十分昂貴，因此利用FFT算法是很有必要的。接下來介紹一種有效的線性迭代法-共軛梯度法（CGS）
三、數值例子
- case $I$(1D) ： 真解：
$$ u = \sin(x),\qquad x\in( 0,\pi ), $$
分別應用直接法和FFT方法的實驗結果見下圖
- case $II$(2D) ： 真解：
$$ u = \sin(x)\sin(y),\qquad (x,y)\in( 0,\pi )^2, $$
分別應用直接法和FFT方法的實驗結果見下圖
從數值實驗結果可以直觀的看出，FFT的計算效率是驚人的！