使用傅里叶变换（包括差分傅里叶和傅里叶谱方法）及切比雪夫谱方法（配点）法求解PDE

傅里叶变换及切比雪夫谱方法（配点）法求解PDE

以下的内容介绍的是傅里叶谱方法求解PDE、切比雪夫谱方法求解PDE（这里指的是配点法）以及一种先进行差分离散，再对离散系统的每个变量使用离散傅里叶级数展开的求解PDE的方法。

因为时间仓促，文中的很多公式，我没能自己手打出来，而用简单地截图来替代，因此也导致了本文在排版上比较混乱，图片大小不一。另外，也因为截图的原因，导致符号不一致（因为截自不同的书籍），符号系统略显混乱，希望读者在阅读时心中自明。

本文内容如下：第一，我想谈谈和傅里叶变换有关的一些东西，以及一种方法使用离散傅里叶变换去求解pde的方法。这和傅里叶谱方法和切比雪夫谱方法是不同的，它是一种离散优先的方法。第二，就是所谓的傅里叶谱方法求解PDE，这和那个“高阶有限元方法”的那个谱方法不尽相同。第三，我想简单地谈一下chebyshev谱方法。

差分离散的傅里叶方法求解PDE

一般的傅里叶变换的定义是：

我们倾向于实用角速度，而不是频率，也就是说：

那么，傅里叶变换就变为了：

傅里叶变换有一些性质，比如说：

如果我用双向箭头来表示傅里叶变换，即：

那么，我们有：

更重要的是，如果我定义卷积（convolution）如下：

那么，有一个卷积定理：

并且，还有Parserval定理：

Parserval定理告诉我们，傅里叶变换在时域（频域）上是平方可积的，那么在（频域）时域上也是平方可积的。当然，傅里叶变换还有很多其他的性质，我不想把它们全都列出来。

下面，我想介绍一个很重要的东西，即香农采样定理。有了它，离散傅里叶变换也就成了顺理成章的事情了。

我们用Δ\DeltaΔ来表示时域上的采样时间间隔，我们能都达到一系列离散的采样点：

我们定义一个特别的频率fcf_cfc，被称为Nyquist critical frequency（奈奎斯特），

它只是和时间间隔Δ\DeltaΔ是有关系的。

什么是采样定理呢？它是说：

我们用ω\omegaω来替代2πf2\pi f2πf，那么写成中文的形式，有：

这个告诉我们什么呢？我们能让采样间隔足够小，以至于f^\hat ff^的支集包含在Nyquist 频率之间，那么采样就没有信息损失。

这个定理非常非常地重要，就像泰勒展开在数值计算中的作用一样。

那么，如果这个条件不被满足，那么会发生什么呢？alising，中文中有的人称为“混叠”，也有叫“截频”之类的，我不知道他们是怎么翻译的。

如下图所示：

那个aliased的那条线其实啊，就是用原来的那个(12.1.3)得到的，只不过不是真正的傅里叶变换之后频域上的取值。其实，误差会以某种方式被推到并且累积到Nyquist区间中（其实有个两者关系的一个表达式）。在我的观点中，Nyquist采样定理是连接连续和离散的桥梁。

下面我们来谈一谈离散傅里叶变换。

我们取：

为了方便，我们假设N是偶数。那么，让我来寻找在傅里叶变换之后，函数值在下述这些点上的估计：

那么，如果我们的采样满足前面的采样条件，且h(t)h(t)h(t)的非零值都集中在一个有限的区间内的话，我们做一个傅里叶变换的离散，能得到以下的公式。可以看到，这个约等号可以使用采样定理h(t)h(t)h(t)的表达式得到，因为那个类似于sinc函数一样的东西，做个傅里叶变换得到的就是指数函数。从另一个角度看，我们通过梯形求积分公式，也能得到类似的表达。

那么，我们再令：

有：

(12.1.7)就是所谓的离散傅里叶变换。容易看到，离散傅里叶变换后的值是周期性的，有周期为N。为了保持一致，我们喜欢变换后的点的下标n从0到N-1。但是从（12.1.5），我们能知道，频率的值是从-N/2 到 N/2。所以呢，有：

DFT有一些性质，和连续傅里叶变换是类似的，没有必要去提它了。

最后，我们定义反（inverse）傅里叶变换如下：

傅里叶变换的用途是什么呢？它可以帮我们求解微分方程，我等等会提到这个。通过以上的讨论，我们知道，对于周期函数，我们可以通过采样，做离散的傅里叶展开，那么对于非周期边界条件呢？

下面，我想介绍一下，快速正弦变换，和快速余弦变换。在如下所示的情况下，我们可以用sine变换，和consine变换。

首先，我需要介绍一下什么是sine变换。

对于fjf_jfj，首先可以关于j=N点，做一个twice的奇延拓，以至于，我们有：

然后，我们再用傅里叶变换：

我们可以把这个式子分成两部分，后面一部分，做一个替换：

把两部分合成一部分，我们有：

当然，我们也能够通过fft算法来进行离散sine变换。

那么离散余弦变换是什么呢？和离散正弦变换的推导相同，此时我们只是在N处做偶延拓，使得：

相同的推导过程，我们能得到：

这就是离散余弦变换。它也有一些快速的算法通过fft。离散余弦变换的另一种形式是：

这里的一撇表示当k=0的时候，那一项前面有个系数0.5。它是原有序列关于点N-1/2做偶延拓得到的。

二维的傅里叶变换定义为：

容易看到，它只是在不同的方向上做一维的傅里叶变换的结果。

到这里，我们已经讨论了离散傅里叶变换，离散正弦变换和离散余弦变换。下面我想说一种使用FFT求解PDE的一种方法，它不同于傅里叶谱方法，和chebyshev谱方法，它是离散优先的一种方法。快速傅里叶方法直接可用，当方程空间表达项前面的系数是常数的时候。

我们以椭圆方程的例子举例。对于问题：

由有限差分方法得到：

重新唤起二维的离散傅里叶变换：

带入到差分问题。我们能得到：

从这里，我们可以看到，我们只需要做一次傅里叶变换，做一次除法，接着傅里叶反变换，我们就完成了pde的求解。

以上的过程对周期性边界是有效的。也就是说：

对于齐次的dirichlet边界变换问题，我们使用sine展开：

和前面相似，过程是：

非齐次边界条件怎么办？书上有一些技巧去处理它。这里不再提了。有兴趣的读者，可以翻翻书。

对于齐次的纽曼边界条件来说，我们使用cos展开：

对于以下的一个问题，我们可以设计一个傅里叶变换解法。我们使用cos变换展开。

时间上，我使用隐式欧拉格式。

离散我们的问题，可以得到：
uj,ls+1=uj.ls+Δth2(uj−1,ls+1+uj+1,ls+1+uj,l−1s+1+uj,l+1s+1−4uj,ls+1)u_{j,l}^{s+1}=u_{j.l}^s+\frac{\Delta t}{h^2}(u_{j-1,l}^{s+1}+u_{j+1,l}^{s+1}+u_{j,l-1}^{s+1}+u_{j,l+1}^{s+1}-4u_{j,l}^{s+1})uj,ls+1=uj.ls+h2Δt(uj−1,ls+1+uj+1,ls+1+uj,l−1s+1+uj,l+1s+1−4uj,ls+1)

替换cos变换进去，cos变换的定义在这。容易知道：
uj−1,ls+1+uj+1,ls+1=2J2L∑m=0J′′∑n=0L′′u^m,ns+1cos(πlnL)[cos(πjmJ)cos(πmJ)+sin(πjmJ)sin(πmJ)+cos(πjmJ)cos(πmJ)−sin(πjmJ)sin(πmJ)]=2J2L∑m=0J′′∑n=0L′′u^m,ns+1cos(πjmJ)cos(πlnL)2cos(πmj)u_{j-1,l}^{s+1}+u_{j+1,l}^{s+1} = \frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi ln}{L}) \\ [cos(\frac{\pi j m}{J})cos(\frac{\pi m}{J})+sin(\frac{\pi j m}{J})sin(\frac{\pi m}{J})\\ +cos(\frac{\pi j m}{J})cos(\frac{\pi m}{J})-sin(\frac{\pi j m}{J})sin(\frac{\pi m}{J})]\\ =\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L})2cos(\frac{\pi m}{j}) uj−1,ls+1+uj+1,ls+1=J2L2m=0∑J′′n=0∑L′′u^m,ns+1cos(Lπln)[cos(Jπjm)cos(Jπm)+sin(Jπjm)sin(Jπm)+cos(Jπjm)cos(Jπm)−sin(Jπjm)sin(Jπm)]=J2L2m=0∑J′′n=0∑L′′u^m,ns+1cos(Jπjm)cos(Lπln)2cos(jπm)

可以得到：
uj−1,ls+1+uj+1,ls+1+uj,l−1s+1+uj,l+1s+1−4uj,ls+1=2J2L∑m=0J′′∑n=0L′′u^m,ns+1cos(πjmJ)cos(πlnL)2[cos(πmj)+cos(πnL)−2]u_{j-1,l}^{s+1}+u_{j+1,l}^{s+1}+u_{j,l-1}^{s+1}+u_{j,l+1}^{s+1}-4u_{j,l}^{s+1}\\ =\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L})2[cos(\frac{\pi m}{j})+cos(\frac{\pi n}{L})-2]uj−1,ls+1+uj+1,ls+1+uj,l−1s+1+uj,l+1s+1−4uj,ls+1=J2L2m=0∑J′′n=0∑L′′u^m,ns+1cos(Jπjm)cos(Lπln)2[cos(jπm)+cos(Lπn)−2]

所以，我们有：

2J2L∑m=0J′′∑n=0L′′u^m,ns+1cos(πjmJ)cos(πlnL)=2J2L∑m=0J′′∑n=0L′′u^m,nscos(πjmJ)cos(πlnL)+Δth2[2J2L∑m=0J′′∑n=0L′′u^m,ns+1cos(πjmJ)cos(πlnL)2[cos(πmj)+cos(πnL)−2]]\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L})\\ =\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L}) \\ +\frac{\Delta t}{h^2}[\frac{2}{J}\frac{2}{L}\sum_{m=0}^J\prime\prime \sum_{n=0}^L\prime\prime \hat u^{s+1}_{m,n}cos(\frac{\pi j m}{J})cos(\frac{\pi ln}{L})2[cos(\frac{\pi m}{j})+cos(\frac{\pi n}{L})-2]]J2L2m=0∑J′′n=0∑L′′u^m,ns+1cos(Jπjm)cos(Lπln)=J2L2m=0∑J′′n=0∑L′′u^m,nscos(Jπjm)cos(Lπln)+h2Δt[J2L2m=0∑J′′n=0∑L′′u^m,ns+1cos(Jπjm)cos(Lπln)2[cos(jπm)+cos(Lπn)−2]]

最后，我们得到（为什么能去掉sum符号？因为e指数表达式一组基）
u^m,ns+1=u^m,ns+Δth2[u^m,ns+12[cos(πmj)+cos(πnL)−2]]\hat u^{s+1}_{m,n} = \hat u^{s}_{m,n} +\frac{\Delta t}{h^2}[\hat u^{s+1}_{m,n}2[cos(\frac{\pi m}{j})+cos(\frac{\pi n}{L})-2]]u^m,ns+1=u^m,ns+h2Δt[u^m,ns+12[cos(jπm)+cos(Lπn)−2]]

也就是说，

u^m,ns+1=u^m,ns/{1−2Δth2[cos(πmj)+cos(πnL)−2]}\hat u^{s+1}_{m,n} = \hat u^{s}_{m,n}/\{1 - \frac{2\Delta t}{h^2}[cos(\frac{\pi m}{j})+cos(\frac{\pi n}{L})-2]\}u^m,ns+1=u^m,ns/{1−h22Δt[cos(jπm)+cos(Lπn)−2]}

思考：这里看起来和sin变换的出来的结果是一样的，这是一个巧合吗？是的。不过，我们这里对uuu做的是cos变换，故而整个结果还是不一样的。

下面，我想通过一个一维的例子说说所谓的傅里叶谱方法求解PDE。

过程是：

那么，我们可以用龙格库塔格式求解了。也就是：

接下来，我想介绍得是chebyshev谱方法。它有点像差分方法，只不过它用的是切比雪夫点，而不是均匀的网格点。

什么是切比雪夫点呢？

我们可以将单位圆上半圆周平均分为N份，等分节点在x轴上的投影就是切比雪夫点。

下面我想介绍的是使用切比雪夫求导矩阵去求一些函数在切比雪夫点上的导数。

想要做的就是：给你一个函数在切比雪夫点上的值，你需要寻找这些点上的导数值。

我们得到离散导数值可以通过两步，即先插值，得到插值多项式，在求导，去这些点上的值：

这个算子是线性的，所以可以写为：

一个矩阵乘上一个向量。先考虑简单的情况。

如果N=1N=1N=1，那么插值点为x0=1,x1=−1x_0=1,x_1=-1x0=1,x1=−1。拉格朗日插值为：

那么求导，得到：

所以呢，求导矩阵为：

考虑N=2N=2N=2的情况，插值节点为1 、0、 -1，拉格朗日插值为：

导数（derivative）为：

这一次，求导矩阵为：

更一般，我们有一个定理。

有了谱求导矩阵，我们就可以使用如下近似（approximation）：

如何使用切比雪夫求导矩阵去求解PDE呢？让我们看一些简单的例子。

考虑一维泊松问题：

通过离散，它可以写为：

如何处理狄利克雷边界条件？这里有N-1个未知数（unknown number），但是有 N+1 个方程。条件太多啦。怎么处理？如下图：

我们先去掉求导矩阵的第一行和最后一行，以及右端项f的第一个值和最后一个值。再由狄利克雷边界条件，删掉求导矩阵的第一列和和最后一列，以及u的第一个值和最后一个值。我们用D~,u~\tilde D,\tilde uD~,u~来表示处理后的结果。最后，求导结果就是：

如果区间不是从-1到1，我们只要在求导矩阵前面乘以一个缩放因子（Scale Factor）。

另一个例子是：

让u~\tilde uu~表示除边界外所有函数的值，即：

D~\tilde DD~也是一样的意思。那么我们有这样的近似：

那么纽曼边界条件怎么办？这也有一个小例子。

边界条件可以写为：

我们可以交换求导顺序（为什么可交换？）：

最后，我们能得到这样一个问题：

我们可以使用Dn2D_n^2Dn2去表示∂2/∂x2\partial ^2/\partial x^2∂2/∂x2。在时间方向上，我们可以用龙格库塔格式。问题是，边界条件怎么办？

not so easy！

如果一个问题有纽曼边界条件，我们有：

在u′(±)=0u'(\pm)=0u′(±)=0的条件下：

到这，我们将边界上的狄利克雷约束条件转化为了简单的边界值表达。

回到原问题，我们将∂u/∂t\partial u / \partial t∂u/∂t视作uuu，那么我们得到∂u/∂t\partial u / \partial t∂u/∂t在边界上的一个表达，直接替换这个表达到切比雪夫离散的右端项上，这样我们时间步（time-stepping）上的格式能够进行。对于切比雪夫离散矩阵，我们将边界上的值替换为其内部（interior）值得一个线性表出。

现在，让我们把这种方法应用到我们的问题（之前写过）上，即：

ut−Δu=0inΩ=[0,1]2∂u∂η∣∂Ω=gη∣∂Ωu(0,x,y)=g(x,y)u_t - \Delta u = 0 \ in \ \Omega=[0,1]^2\\ \frac{\partial u}{\partial \eta}\mid _{\partial \Omega}=\frac{g}{\eta}\mid _{\partial \Omega}\\ u(0,x,y) = g(x,y)ut−Δu=0 in Ω=[0,1]2∂η∂u∣∂Ω=ηg∣∂Ωu(0,x,y)=g(x,y)

容易convert这个问题为：

ut−Δu=0inΩ=[0,1]2∂∂η∂u∂t=0u(0,x,y)=g(x,y)u_t - \Delta u = 0 \ in \ \Omega=[0,1]^2\\ \frac{\partial }{\partial \eta}\frac{\partial u}{\partial t}=0\\ u(0,x,y) = g(x,y)ut−Δu=0 in Ω=[0,1]2∂η∂∂t∂u=0u(0,x,y)=g(x,y)

然后，我们在不同的方向上，对如上的问题如法炮制，我们能够很快地求解这个问题。

下一文章，我想介绍的是，傅里叶级数（三角形式和指数形式）、连续傅里叶变换、离散时间傅里叶变换、傅里叶变换的所需要的条件（对函数的要求），以及他们之间的一个关系。敬请期待。