# 傅里叶谱方法及其工程实现

## 傅里叶谱方法求解PDE

#### 快速傅里叶变换

u ^ ( k ) = ∫ − ∞ + ∞ u ( x ) e − i k x d x u ( k ) = 1 2 π ∫ − ∞ + ∞ u ^ ( x ) e i k x d x \hat u(k)=\int_{-\infty}^{+\infty}u(x)e^{-ikx}dx\\ u(k)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\hat u(x)e^{ikx}dx

#### DFT和FFT

DFT，离散傅里叶变换，就是傅里叶变换的离散化。简单地说，就是把一个N长的向量，通过乘以一个MxN的矩阵，变成M长的向量（一般可以去M=N）。重要的是，这个变换矩阵是什么呢？其实它是一个范德蒙德矩阵，所以DFT其实就可以写为：

u ^ k = ∑ j = 0 N − 1 e i 2 π k M ⋅ j ⋅ u j , j = 1 , 2 , . . . , N \hat u_k = \sum\limits_{j=0}^{N-1}{e^{i\frac{2\pi k}{M}\cdot j}\cdot u_j},j=1,2,...,N

#### 傅里叶谱方法

F [ u ( n ) ( x ) ] = ( i k ) n F [ u ( x ) ] {F[u^{(n)}(x)]}=(ik)^nF[u(x)]

F − 1 { F [ u ( n ) ( x ) ] / ( i k ) n } = u ( x ) F^{-1}\{F[u^{(n)}(x)]/(ik)^n\}=u(x)

#### 一个简单的数值算例

( Δ u ) ^ ( k ) = ∫ Ω Δ u ( x ) e − i k x d x = ∫ ∂ Ω ∂ u ∂ η e − i k x d s + i ∣ k ∣ ∫ Ω ∇ u ( x ) e − i k x d x = i ∣ k ∣ ∫ Ω ∇ u ( x ) e − i k x d x = i ∣ k ∣ ∫ ∂ Ω u ( x ) e − i k x d s + ( i ∣ k ∣ ) 2 ∫ Ω u ( x ) e − i k x d x = − ∣ k ∣ 2 ∫ Ω u ( x ) e − i k x d x = − ∣ k ∣ 2 u ^ ( x ) \hat{(\Delta u)}(k)=\int _{\Omega}\Delta u(x)e^{-ikx}dx\\ =\int_{\partial \Omega}\frac{\partial u}{\partial \eta}e^{-ikx}ds+i|k|\int_{\Omega }\nabla u(x)e^{-ikx}dx\\ =i|k|\int_{\Omega }\nabla u(x)e^{-ikx}dx\\ =i|k|\int_{\partial \Omega}u(x)e^{-ikx}ds+(i|k|)^2\int_{\Omega }u(x)e^{-ikx}dx\\ =-|k|^2\int_{\Omega }u(x)e^{-ikx}dx= -|k|^2\hat u(x)

Δ u ^ ( k ) = − ( k 1 2 + k 2 2 ) u ^ ( x ) \hat{\Delta u}(k) = -(k_1^2+k_2^2)\hat u(x)

∂ u ^ ∂ t + ( k 1 2 + k 2 2 ) u ^ = 0 u ^ ( 0 ) = g ^ \frac{\partial \hat u}{\partial t}+(k_1^2+k_2^2)\hat u = 0\\ \hat u(0) = \hat g

function T = Runge_Kutta_fft_spectral(D,h,tau)
%谱方法求解二维热传导方程
L=1; N=round(1/h);%划分
x=L/N*[-N/2:N/2-1]; y=x;%拆成两个维度来
kx=(2*pi/L)*[0:N/2-1 -N/2:-1]; ky=kx;%对应的频域
[X,Y]=meshgrid(x,y);
[kX,kY]=meshgrid(kx,ky);
K2=kX.^2+kY.^2;%非顺序下的K方
%初始条件
D0 = D(1:end-1,1:end-1);
u0 = D0;
ut0 = fft2(u0);
ut= ut0; ut_flatten = ut(:);%flatten
K2_flatten = K2(:);
%求解
t=[0 tau];
%options = odeset('RelTol',1e-5,'AbsTol',1e-8);
[t,utsol]=ode45(@(t,ut_flatten) heat_equation_2D(t,ut_flatten,K2_flatten),t,ut_flatten);
T_temp = ifft2(reshape(utsol(end,1:N^2),N,N));

T_temp = [T_temp T_temp(:,1)];
T = [T_temp;T_temp(1,:)];%周期性边界扩充
end
function dut=heat_equation_2D(t,ut_flatten,K2_flatten)
dut=[-K2_flatten.*ut_flatten];
end


a = 1 a=1 ，初始条件为：

#### 其他一些东西

• 滤波法
对于刚性比较大的微分方程，所以的刚性大，就是扰动太大，变化太快，我们可以使用滤波法来降低刚性。

• 谱求导矩阵
谱方法还可以用来做插值。使用谱求导矩阵，可以确定插值过程中的插值函数，并将其用于计算离散数据的各阶导数、偏微分方程（组）的数值求解。

• 切比雪夫谱方法
不管是傅里叶谱方法还是谱求导周期插值方法，在进行数值求导时都暗含了周期性边界条件，很有局限性。对于一般的第一类、第二类和第三类边界条件问题就不好处理。所以就有了切比雪夫谱方法。也就是说，切比雪夫谱方法可以处理非周期区间内的问题。

01-27 2万+

06-15 6万+
12-06 1815
09-09
09-11
11-29 3万+
01-09 10万+