Matlab谱方法：傅里叶谱解法的步骤

待求解的偏微分方程的普遍形式为：

其中，u(x,t)为x、t的函数，L代表线性算符（linear operator），N(u)为非线性项（nonlinear terms）。已知初始条件为u(x,t₀)，在周期性边界条件下求某一时刻的u(x,t)。对于这样的问题，总是可以通过变量代换的办法将等号左边对t的n阶导数降为1阶，详见1.1.6小节。所以下面以等号左边是∂u/∂t的偏微分方程为例进行讨论，这与式（3-11）降阶后的偏微分方程组求解方法大同小异。

在x域上对以下方程做傅里叶变换：

得到：

其中，û(k,t)代表u(x,t)在x域上的傅里叶变换。为了展现对线性算符L及非线性项N(u)做傅里叶变换的细节，这里以L=a·∂²/∂x²+b·∂/∂x+c、N(u)=u³+u³·∂²u/∂x²+f(x)·∂u/∂x为例进行分析，a、b和c分别为常数。

对线性部分，有F[Lu]=[a(ik)2+b(ik)+c]û，则式（3-13）中的α(k)=a(ik)2+b(ik)+c=­ak²+ibk+c。对非线性部分，需要将u、∂u/∂x和∂²u/∂x²写为F^-1[û]、F^-1[ikû]和F^-1[(ik)²û]，则有F[N(u)]=F[u³+u³·∂²u/∂x²+f(x)·∂u/∂x]=F{F^-1[û]³+F^-1[û]³·F^-1[(ik)²û]+f(x)·F^-1[ikû]}。

这样，式（3-12）中u(x,t)对x的偏导数就都通过傅里叶变换及逆变换简化为式（3-13）中û(k,t)和k的代数运算了，然后再将û和k离散化，偏微分方程就简化成了常微分方程组。数值计算最直接的方法就是调用Matlab的ode系列函数，优先选择ode45函数，最后用ifft函数将频域上的计算结果变换回待求的u(x,t)。

图3-3 一维情况下的傅里叶谱方法计算过程(www.xing528.com)

综上，利用傅里叶谱方法数值计算偏微分方程（组）步骤如图3-3所示，具体过程如下：

（1）对于形如式（3-11）的偏微分方程（组），通过变量代换将其中的∂ⁿu/∂tⁿ降为1阶导数，若n=1，略过此步。

（2）在x域上对偏微分方程（组）做傅里叶变换，则线性项中的∂ⁿu/∂xⁿ直接变为û，并利用、代换非线性项中的u、∂ⁿu/∂xⁿ，得到关于的微分方程（组）。

（3）用时间步进法（如龙格-库塔法等）数值计算离散化的关于的微分方程组，默认边界条件为周期性边界条件，初始条件。

（4）将上一步得到的结果从频域变换回空域，即：。离散傅里叶变换及逆变换由fft、ifft函数实现。

若方程在空间上再增加一个维度，式（3-12）的等号右边包含u(x,y,t)的∂ⁿ/∂xⁿ和∂ⁿ/∂yⁿ，则傅里叶谱方法的计算步骤基本不变，如示意图3-4所示。唯一不同的是需要对方程做二维傅里叶变换，使用fft2、ifft2函数代替fft、ifft函数，并利用∂ⁿ/∂xⁿ→(ik_x)n、∂ⁿ/∂yⁿ→(ik_y)n和关系式（3-14）进行转化，其中k_x、k_y为x、y方向的波数。式中F[]、F^-1[]代表二维傅里叶变换及逆变换，u上加“~”和“^”代表先后对函数u做x、y方向上的傅里叶变换（即二维傅里叶变换）。

图3-4 二维情况下的傅里叶谱方法计算过程

为了更贴近实际地说明具体计算细节，后面章节将给出形如式（3-11）和（3-12）的偏微分方程（组）的求解实例。

需要强调的是，使用傅里叶谱方法求解偏微分方程（组）隐含着周期性边界条件。以序列u₁,…,u_j,…,u_N为例，在周期性边界条件下，可以等效认为有u_mN+j=u_j的关系（m为任意整数），也就是说一端边界处的函数值将对另一端边界处的函数值产生影响。周期性边界条件可以将具有时空周期性的物理问题简化为单元进行处理，但对于一些特定的非周期性问题，则需要修改其他条件（计算区间的范围、参数、初始条件等）来确保边界处的函数值恒为0或某一特定常数，以排除相邻周期间的干扰，得到正确结果。