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

傅里叶谱方法求解PDE

快速傅里叶变换

用一维举例，在$(-\infty ,+\infty )$有定义且绝对可积，它的傅里叶变换及其逆变换为：
$$\begin{gathered} \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 \end{gathered}$$
傅里叶变换及其逆变换称为傅里叶变换对，傅里叶变换对不是唯一的。两个定义式子中的系数可以随意修改，只要他们的乘积是$\frac{1}{2\pi }$就行了。当然定义式中的$e^{-ikx}$和$e^{ikx}$也是可以互换的，它们叫做记分变换的核。
傅里叶变换，大家喜欢说是从时域（空域）到频域上的变换。看图稍微理解一下，不再细述。因为不是关键。

频域的基本单元也可以理解为一个始终在旋转的圆。简单地说就是，如果傅里叶变换表示的是从时域到频域上的变换，那么$k$表示的就是频率（角速度，有人用$\omega$来表示），$\hat{u}(x)$表示的就是频谱。类似，如果$x$表示的是空间坐标，那么傅里叶变换就是从空域到频域上的一个变换。那么$k$表示的就是波数（$2\pi$ 长度上的波长的个数，代表空间上的变化速度，暗含频率的意思），此时$\hat{u}(k)$称为波数谱。

简单地说，傅里叶变换及其逆变换就是从时域（空域）到频域之间相互转换的工具。

DFT和FFT

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

那么这个${\omega }_0,{\omega }_1...{\omega }_{M-1}$是个什么玩意儿呢？它其实就是复单位圆周上的M个等分点位置。如下：

一言以蔽之，要将一个N长的向量变成M长的向量，将复平面的上单元圆分成M分，取第1个等分点，乘方出来N个数，和原来的N长向量做內积，得到傅里叶变换出来的第1个数，依次类推，每个等分点都能出来一个数，那么M个等分点就能出来M个数。

我想到这，我应该把离散傅里叶变换将清楚了。非要用公式来表达的话，可以写成这样。

$${\hat{u}}_k=\sum _{j=0}^{N-1}{e}^{i\frac{2\pi k}{M}\cdot j}\cdot u_j,j=1,2,...,N$$

上面提到的都是一维的傅里叶变换，二维傅里叶变换其实就是在两个维度上分别做傅里叶变换。

在实际计算中，人们不直接采用上述离散傅里叶变换的定义式进行计算，而是使用快速傅里叶变换（Fast Fourier Transform），我们可以简单地先把快速傅里叶变换认为是一种做离散傅里叶变换的快速算法。

快速傅里叶变换被誉为二十世纪最伟大的算法之一。DFT算法的运算量为$O(N^2)$，而快速傅里叶变换将运算量降为了$O(NlogN)$。但是，为了获得较高的运算速度，待变换序列的元素数量必须是$2^n$形式个。

傅里叶谱方法

傅里叶谱方法的思路：对PDE方程进行FFT变换，得到只对时间微分的常微分方程组。然后，使用一个基于时间步的数值格式来解ODE方程组。

在matlab中，采用的离散傅里叶变换的定义为（下面只考虑N=M）：

因为matlab的下标是从1开始，这样定义式比较好的。注意到上述定义中的归一化系数也可以其他选择，但它们的成绩必须为$\frac{1}{N}$。并且，由此定义，我们可以看到${\hat{u}}_k={\hat{u}}_{k+N},u_j=u_{j+N}$，也就是说，离散傅里叶变换，隐含着周期性边界条件。

在matlab中，调用傅里叶变换和逆变换的函数为$\text{fft(u)}$和$\text{ifft(u)}$。若u为以矩阵，则fft返回矩阵每一列的离散傅里叶变换。若要计算矩阵u的每一行的离散傅里叶变换，只需调用fft(u,[],2)即可，它比先转置做傅里叶变换，再转置回来要快一些。

需要极其注意的一点是，fft作用于离散函数向量值$u$，它所对应的$\hat{u}$，对应的频域上的定义域（频率$k$）不是从小到大排列的，而是在0位置砍断之后的一个左右平移交换。什么意思呢？下面需要仔细说说这个东西，这个对傅里叶谱方法的实现尤为重要。

空域（时域）上的序列$u_1,u_2,\dots ,u_N$在x轴上的间隔是$\frac{L}{N}$，相应的横坐标为：

为了方便，这里的区间我们取为$[-L/2,L/2]$，因为我们分成N分，产生了N+1个点，由于周期性边界的假设，我们抛弃最右边那个点，得到了N个点。那么$u$序列做傅里叶变换的输出序列${\hat{u}}_1,{\hat{u}}_2,\dots ,{\hat{u}}_N$，其对应的频率的顺序排列为：

注意到，这正是在$2\pi N/L$的区间长度上的等分点，在中间割开之后，左右进行了互换。

那么在变换过去的频域上的区间长度和两点间隔是多少呢？看下面这个表格。

怎么记？有这样一个规律，上表格的对角线相乘是横等于$2\pi$的，所以只要知道了空域坐标的分割情况，就知道频域坐标的分割。为什么这样？现在我们在乎的工程实现，其中的道理先不用深究。

在matlab当中啊，有fftshift函数和ifftshift函数来调整输出的$\hat{u}$的顺序，使其与顺序的频域坐标对应上。比如fftshift(fft(u))，就使得，输出的$\hat{u}$对应的频域坐标是按从小到大（递增）的顺序排列的。当然，在程序中我们考虑变频域坐标是更为方便的，否则在做逆变换的时候，就需要将$\hat{u}$变过去，再变回来。什么意思呢？后面在程序中就能看到这一点。ifftshift函数用于进行与fftshift函数相反的操作。

二维的傅里叶变换matlab调用语法为fft2(u)，由于二维傅里叶变换是做两次一维傅里叶变换，所以其等价于fft(fft(u),[],2)，即先对列做傅里叶变换，再对行做傅里叶变换。同样，二维傅里叶变换，变换过去的值，对应的频率（空间坐标）也不是对应的。那么，它也有一个移位函数fftshift(u)。易知，这个函数是将$\hat{u}(x,y)$关于x轴左右互换，再关于y轴上下互换得到的，其实就是将其一、三象限交换，二、四象限交换并返回。其实就是先左右倒一下，再上下倒一下。

也就是说，做傅里叶变换之后，低频部分到了两端，高频部分到了中间。在实际中，我们常常使用abs函数求绝对值，来得到频谱的振幅。

下面的一个问题是，求解PDE的时候为什么要做傅里叶变换？它到底有什么用呢？

若用F表示傅里叶变换，若满足：在$|x|\to \infty$时，有$u(x)\to 0$，那么由分部积分公式容易得到：

$$F[u^{(n)}(x)]=(ik{)}^{n}F[u(x)]$$

这里的k是频域坐标。导数的傅里叶变换，相当于傅里叶变换穿进去，出来一个$ik$的n方。也就是说，函数的求导运算在傅里叶变换的作用下，可以转化为相对简单的代数运算。也就是要求n阶到，只需要做个傅里叶变换，乘$(ik{)}^n$，在变回去就完事了。这是其在求上的意义。实际操作当中，我们要注意到这个$k$和$F[u(x)]$不是直接对应的，差了一个移位关系，所以在做乘法的时候要注意。对于高维，比如说二维，这个表达式也是类似的，我们可以直接根据定义，对两个方向分别做傅里叶变换和分布积分得到结果，也可以直接把$x$视为一个向量，利用格林公式来处理，得到的是相同的结果。

同样地，对上式进行稍加处理，得到：

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

那么傅里叶变换就具有了求不定积分的意义。用这个来计算积分，在实际操作中，可能会出现$k=0$，而导致式子出现无穷大，通常可以通过将$k$修改为一个很小的数，或者分母加上很小的数来解决。

傅里叶变换还可以用来求解PDE，什么意思呢？我们知道非定常的发展PDE，方程中既带有对时间的导数，又带有对空间的导数，从上面可以看到，对于空间的导数，我们关于空间变量做一个傅里叶变换，就不含对空间（这里的空间指的是频率空间）的导数了，而对时间的导数项，做傅里叶变换，不影响空间坐标上的变化。所以，通过傅里叶变换，我们可以将$u,\partial u/\partial x,{\partial }^2/\partial x^2\dots$变为$\hat{u},ik\hat{u},(ik{)}^2\hat{u}$。容易知道，通过这样的变化，PDE就变成了一个只关于时间变量t的ODE。然后，我们再用一般的ODE的数值方法求解即可，比如用matlab内置ode45，一般使用时ode相关函数时，优先选择这个。最后，把结果，再做一个傅里叶逆变换，变回来，就完事了。是不是很简单呢？

需要注意的一点是，傅里叶谱方法求解PDE一般默认边界条件为周期性的边界添加，也就是说一端边界处的函数值将对另一端边界处的函数值产生影响。对于一些非周期性的文艺，我们如何确保边界处的值对整体影响不大或者能确保边界处的值一直为一个特定的常数，是一个值得思考的问题。

一个简单的数值算例

一个简单的数值算例如下：

这里的$R_0=1/4$，我们计算$t=3/256$时刻的$u$值。假设它是周期性边界条件，将$\Omega$扩充到$R^2$上，那么有：

$$\begin{gathered} \widehat{(\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) \end{gathered}$$

注意到，这里有一个问题，就是要使上式成立，要保证在$x\to -\infty$时，有$u(x)\to 0$，并且趋于0的速度要大于指数速率，才能保证边界项为0。

如果对二维格林公式的操作不熟，我们也可以将$u$写为$u(t,x,y)$，然后对于$\Delta u(t,x,y)$分别对$x$和$y$做两次一维的傅里叶变换，结果是一样的。即：
$$\widehat{\Delta u}(k)=-(k_1^2+k_2^2)\hat{u}(x)$$

那么，记$u(0,x,y)=g(x)$对上述问题，做傅里叶变换，就得到了：

$$\begin{gathered} \frac{\partial \hat{u}}{\partial t}+(k_1^2+k_2^2)\hat{u}=0 \\ \hat{u}(0)=\hat{g} \end{gathered}$$
对于这样一个ODE问题，我们可以设计欧拉格式、龙哥库塔等方法求解。写出来一个matlab脚本如下：

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

程序比较好懂，注意看到的是我们的$k_x,k_y$的取值为(2*pi/L)*[0:N/2-1 -N/2:-1]，是为了计算$(k_1^2+k_2^2)\hat{u}$的时候，频域坐标对应的频谱值是正确的。

以上就是求解PDE的傅里叶谱方法，还有一个需要注意的地方是，但时间上的导数不只一阶的时候，我们往往通过引入一些新的变量（正如处理高阶ODE那样），将高阶的ODE变成一阶的ODE方程组。在做傅里叶变换，从而使用傅里叶谱方法。举个简单的例子。

对于无界区域上的二维波动方程：

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

引入函数$v=\frac{\partial u}{\partial t}$，进行降阶，得到：

对上式等号两边做傅里叶变换，得到：

再用数值方法求解此ODE方程组即可。

其他一些东西

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

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

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

谱方法

谱方法和有限元法的思想很类似，差别在于：谱方法以一系列全局连续的函数（可以是三角函数，也可以是多项式）的叠加来近似真实解，而有限元法则是使用分片的简单函数的叠加来近似近似真实解。即，有限元的插值函数，只在该单元内作用。

谱方法方法是解偏微分方程的一种数值方法。其要点是把解近似地展开成光滑函数(一般是正交多项式)的有限级数展开式，即所谓解的近似谱展开式，再根据此展开式和原方程，求出展开式系数的方程组。对于非定常问题，方程组还同时间有关。

谱方法实质上是标准的分离变量技术的一种推广。一般多取切比雪夫多项式和勒让德多项式作为近似展开式的基函数。对于周期性边界条件，用傅里叶级数和调和级数比较方便。谱方法的精度，直接取决于级数展开式的项数。

一般说，谱方法远比普通一、二阶差分法准确。由于快速傅里叶变换之类的技术不断发展，谱方法的运算量越来越少，一般是很合算的。特别是对于二维以上的问题，用差分法计算必须设置足够多的网格点，造成计算量的增加，而用谱方法一般不需取太多的项就可得到较高精度的解。因此谱方法在计算流体力学复杂流场的问题中有广泛应用。