数值方法求解微分方程的研究——基于切比雪夫多项式的谱方法

Research on Solving Differential Equations by Numerical Method—Spectral Methods Based on Chebyshev Polynomials

摘要: 谱方法是处理微分方程的常用方法，本文以理论完善的谱方法为基础，详细介绍了切比雪夫多项式通过S-L问题的由来与切比雪夫多项式的部分性质，并利用这些性质将这些正交多项式作为基对函数进行展开，从而数值求解偏微分方程，我们利用案例来展现其具体的运算过程并验证其方法的有效性。

Abstract: Based on Spectral Method, we introduce the origin of Chebyshev polynomials via S-L problems and some properties of Chebyshev polynomials. With these properties, we use these orthogonal poly-nomials as basic functions to solve partial differential equations. Also, we use examples to show the detailed operations and verify its effectiveness.

1. 引言

偏微分方程是反映未知变量关于时间的导数和关于空间变量的导数之间制约关系的等式。许多领域中的数学模型都可以利用偏微分方程来描述，甚至很多数学、物理和力学分支中的基本方程本身就是微分方程。经过许多年的发展，偏微分方程已经成为当代数学中的一个重要组成部分，是纯数学学科和其他工程学科、自然科学之间的一座桥梁 [1] 。

对于十八、十九世纪建立起来的众多微分方程，数学家们求其显式解往往失败，这种情况促使他们转而证明解的存在性。而研究物理问题过程出现的具体微分方程问题时，往往会使用不同的方法，例如分离系数法、分离变数法，也有许多问题不能严格解出，只能用近似的方法求出满足实际需要的近似程度的解，常用的方法例如变分法、有限差分法。谱方法也很早就出现，它是利用一组基函数将微分方程的解表示出来，但由于计算量大一直没有被广泛使用，直到快速Fourier变换的出现，给谱方法带来了生机，近几十年快速发展，成为有限差分法、有限元法后的第三种数值求解的基本方法。

谱方法最大的魅力就是它具有“无穷阶”收敛性，如果原方程的解无穷光滑，那么用适用的谱方法所求得的近似解将以 $N^{-1}$ 的任意幂次速度收敛于精确解，这里的N为所选取的基函数的个数。这一优点是有限差分法和有限元法无法比拟的，众多的实际应用和数值实例也证实了这一方法的有效性。因此谱方法日益受到重视，但是和差分法与有限元法的研究相比，仍存在一定差距，谱方法在非线性情形下，或问题不是周期情况下，特别是在求解偏微分方程数值解方面也需要更多深入的研究。谱方法理论的日益完善，也被用于处理物理、数学、工程中各种数值问题的处理，比如麦克斯韦方程，也以谱方法为基础，衍生不少新的方法 [2] 。

2. 切比雪夫多项式

2.1. Sturm-Liouville (S-L)问题

微分方程的谱近似都可以视为S-L问题特诊函数的有限展开，方程如下 [3] [4] ：

$-{\left(p{u}^{\prime }\right)}^{\prime }+qu=\lambda \omega uu\in \left(-1,1\right)$

其中p，q，ω都是给定的，它们是实值函数且满足：p是连续可微的，在(−1, 1)区间内是正值且在边界 $x=\pm 1$ 处连续；q在(−1, 1)内连续非负；权重函数ω在(−1, 1)上连续、非负且可积。方程的一些非平凡的解λ称为方程的特征值。方程的解为 $\left({\lambda }_0,u_0\left(x\right)\right),\left({\lambda }_1,u_1\left(x\right)\right),\left({\lambda }_3,u_3\left(x\right)\right),\cdots$ 。

在奇异S-L问题中，代数多项式可以用数值方法进行评估和微分，方程组在(−1, 1)上的解满足正交性，这个性质使得方程的解：如勒让德多项式、切比雪夫多项式，由其不同性质可以用于问题的数值求解。

2.2. 正交多项式

考虑一族定义在(−1, 1)上的多项式 $p_0\left(x\right),p_1\left(x\right),p_2\left(x\right),\cdots$ 若满足

${\left(p_m\left(x\right),p_n\left(x\right)\right)}_{\omega }={\displaystyle {\int }_{-1}^1{p}_m\left(x\right)p_n\left(x\right)\omega \left(x\right)\text{d}x}=0,对一切m\ne n$

那么它是(−1, 1)上的正交多项式系。当权重函数 $\omega \equiv 1$ 时便是最常见的形式

$\left(p_m\left(x\right),p_n\left(x\right)\right)={\displaystyle {\int }_{-1}^1{p}_m\left(x\right)p_n\left(x\right)\text{d}x}=0,对一切m\ne n$

正常我们只在有限维内(即多项式的个数小于等于N)考虑，满足上面条件即在(−1, 1)上两两正交。

2.3. 由S-L问题引出切比雪夫多项式

${\left(\sqrt{1-x^2}{T^{\prime }}_k\left(x\right)\right)}^{\prime }+\frac{k^2}{\sqrt{1-x^2}}{T}_k\left(x\right)=0$

问题令S-L问题中

$p\left(x\right)=\sqrt{1-x^2},q\left(x\right)=0,\omega \left(x\right)=\frac{1}{\sqrt{1-x^2}}$

可以得到

$T_k\left(x\right)=\cos kx,z=\arccos x$

或者写成

$T_k=\cos \left(k{\cos }^{-1}x\right)$

因切比雪夫多项式是由S-L问题推导出的特殊情况，故其内积也满足正交性：即当 $k\ne l$ 时 ${\left(T_m,T_l\right)}_{\omega }=0$ ，且当 $k=l=0$ 时 ${\left(T_m,T_l\right)}_{\omega }=\pi$ ， $k=l\ne 0$ 时 ${\left(T_m,T_l\right)}_{\omega }=\pi /2$ 。写出其具体表达式：

$T_0=1,T_1=\cos z=x,T_2=\cos 2z=2{\cos }^{2}z-1=2x^2-1$

而由等式

$\cos \left(k+1\right)z+\cos \left(k-1\right)z=2\cos z\cos kz$

将切比雪夫多项式带入有

$T_{k+1}-2x{T}_k+T_{k-1}=0$

利用此公式，可以得到

$T_3=2\times x{T}_2-T_1=4x^3-3x$

如果需要，也不难推导出其后面每项的具体展开表达式。

2.4. 切比雪夫展开

类似傅里叶展开，我们也可以将函数用切比雪夫多项式展开成有限个多项式的和的形式，即：

$u\left(x\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{\hat{u}}_k{T}_k(\; x\; )}$

其中

${\hat{u}}_k=\frac{2}{\pi C_k}{\displaystyle {\int }_{-1}^{1}u\left(x\right)T_k\left(x\right)\omega \left(x\right)\text{d}x}$

$C_k=\{\begin{array}{l}2,k=0\\ 1,k\ge 1\end{array}$

2.5. 方程的逼近方法与高斯积分

将有限项：其中 $T_k\left(x\right)$ 为切比雪夫多项式，选其作为试函数(trial function)。从方程近似的角度来看，典型的谱方法有两种：分为在物理空间离散求解的Collocation法与在谱空间进行离散求解的Galerkin法。具体做法Collocation法是在定义域内取点，取的点的数量取决于试函数中待定系数的个数，使得残差在这些点上等于零。取点的位置可以是随意的，但一般均匀取点或者根据分析调整。Galerkin法为通过对基函数在定义域内的加权积分满足原方程，得到一组易于求解的线性代数方程，边界条件自然能够满足。Collocation法适用于非线性问题，Galerkin法适用于线性问题。

在用Collocation法计算时，可以选用点如下搭配点

$x_j=\cos \frac{\pi j}{N}$

选取其不仅可以保证的计算的精确度，而且可以利用快速傅里叶变换使计算更加简便。

在计算展开系数时，通常可以用数值方法高斯积分来简化计算过程。注意下面这个方程组：

$\underset{j=0}{\overset{N}{{\displaystyle \sum }}}{\left(x_j\right)}^k{\omega }_j={\displaystyle {\int }_{-1}^1{x}^k\omega \left(x\right)\text{d}x}$

其中 $x_0<x_1<\cdots <x_N$ 为N + 1阶正交多项式 $P_{N+1}$ 的根， ${\stackrel{\dot{}}{\omega }}_0,\cdots ,{\omega }_N$ 为上线性方程的解，但是 $x_i$ 与 ${\omega }_i$ 不容易找到，便引出了高斯–拉道(Gauss-Radau)、高斯–洛巴托(Guass-Lobatto)积分。而我们选取的搭配点来自高斯–洛巴托积分：对于 $q\left(x\right)=P_{N+1}\left(x\right)+a{P}_N\left(x\right)+b{P}_{N-1}\left(x\right)$ ，选取a，b使得 $q\left(-1\right)=q\left(1\right)=0$ ， $-1=x_0<x_1<\cdots <x_N$ 为 $q\left(x\right)$ 的根，对

$\underset{j=0}{\overset{N}{{\displaystyle \sum }}}{\left(x_j\right)}^k{\omega }_j={\displaystyle {\int }_{-1}^1{x}^k\omega \left(x\right)\text{d}x},0\le k\le N$

有

$\underset{j=0}{\overset{N}{{\displaystyle \sum }}}P\left(x_j\right){\omega }_j={\displaystyle {\int }_{-1}^{1}P\left(x\right)\omega \left(x\right)\text{d}x},任意P\in {\mathbb{P}}_{2N-1}$

2.6. 系数的推导

取前N + 1项对等式两边同时乘以 $T_m\left(x\right)$ 并取内积

${\left(T_m\left(x\right),u\left(x\right)\right)}_{\omega }={\left(\underset{k=0}{\overset{N}{{\displaystyle \sum }}}{\hat{u}}_k{T}_k\left(x\right),T_m\left(x\right)\right)}_{\omega }$

由切比雪夫多项式的正交性，右侧当 $k\ne m$ 时恒为0，k = m时有

${\left({\hat{u}}_m{T}_m\left(x\right),T_m\left(x\right)\right)}_{\omega }={\hat{u}}_m{\left(T_m\left(x\right),T_m\left(x\right)\right)}_{\omega }=\frac{\pi }{2}{c}_m{\hat{u}}_m$

其中

$c_m=\{\begin{array}{l}2,m=0\\ 1,m\ge 1\end{array}$

左侧应用上述高斯积分公式

${\displaystyle {\int }_{-1}^1{T}_m\left(x\right)u\left(x\right)\omega \left(x\right)\text{d}x}=\frac{\pi }{N}\underset{i=0}{\overset{N}{{\displaystyle \sum }}}\frac{T_m\left(x_i\right)u\left(x\right)}{{\bar{c}}_i}$

其中

${\bar{c}}_i=\{\begin{array}{l}2i=0/k=N\\ 11\le i\le N-1\end{array}$

两侧联立可以得到上述展开多项式的系数 [5] [6]

${\hat{u}}_K=\frac{2}{{\bar{c}}_{k}N}\underset{i=0}{\overset{N}{{\displaystyle \sum }}}\frac{1}{{\bar{C}}_i}{u}_i\cos \frac{k\pi i}{N},k=0,\cdots ,N$

2.7. 微分

考虑其展开的前N + 1项

$u_N\left(x\right)=\underset{k=0}{\overset{N}{{\displaystyle \sum }}}{\hat{u}}_k{T}_k(\; x\; )$

则其一阶微分写成

${U^{\prime }}_N\left(x\right)=\underset{k=0}{\overset{N}{{\displaystyle \sum }}}{\hat{u}}_k{T^{\prime }}_k\left(x\right)=\underset{k=0}{\overset{N}{{\displaystyle \sum }}}{u}_k^{\left(1\right)}{T}_k(\; x\; )$

${T^{\prime }}_k=\frac{\text{d}}{\text{d}z}\left(\cos kz\right)\frac{\text{d}z}{\text{d}x}=k\frac{\sin kz}{\sin z}$

同时又有 $2\sin \theta \cos k\theta =\sin \left(k+1\right)\theta -\sin \left(k-1\right)\theta$ ，两边同除 $\sin \theta$

$\Rightarrow 2T_k\left(x\right)=\frac{1}{k+1}{T^{\prime }}_{k+1}\left(x\right)-\frac{1}{k-1}{T^{\prime }}_{k-1}(\; x\; )$

将其带入一阶微分等式，得到

$2\underset{k=0}{\overset{N}{{\displaystyle \sum }}}{\hat{u}}_k{T^{\prime }}_k\left(x\right)=\underset{k=0}{\overset{N}{{\displaystyle \sum }}}{u}_k^{\left(1\right)}\left(\frac{1}{k+1}{T^{\prime }}_{k+1}-\frac{1}{k-1}{T^{\prime }}_{k-1}\right)$

使 ${T^{\prime }}_k$ 对应的系数相等，左侧为 $2{\hat{u}}_k$ ，右侧分别为

$u_{k-1}^{\left(1\right)}\frac{1}{k}$ 与 $-u_{k+1}^{\left(1\right)}\frac{1}{k}$

考虑上式要满足k ≥ 2则有

$2k{\hat{u}}_k=u_{k-1}^{\left(1\right)}-u_{k+1}^{\left(1\right)}\left(k\ge 2\right)$

当k = 1时单独考虑，综合上式总结得到：

$2k{\hat{u}}_k=C_{k-1}{u}_{k-1}^{\left(1\right)}-u_{k+1}^{(\; 1\; )}$

当 $k\ge N$ 时， ${\hat{u}}_k^{\left(1\right)}=0$

当 $0\le k\le N-1$ 时， $C_k{\hat{u}}_k^{\left(1\right)}={\hat{u}}_{k+2}^{\left(1\right)}+2\left(k+1\right){\hat{u}}_{k+1}$

其具体形式为 [7] [8] [9] ：

$C_{N-1}{\hat{u}}_{N-1}^{\left(1\right)}=2N{\hat{u}}_N$

$C_{N-2}{\hat{u}}_{N-2}^{\left(1\right)}=2\left(N-1\right){\hat{u}}_{N-1}$

$C_{N-3}{\hat{u}}_{N-3}={\hat{u}}_{N-1}+2\left(N-2\right){\hat{u}}_{N-2}=\frac{2}{C_{N-1}}N{\hat{u}}_N+2\left(N-2\right){\hat{u}}_{N-2}$

3. 计算案例

3.1. 热传导方程概括

$\frac{\partial u}{\partial t}-\frac{{\partial }^{2}u}{\partial x^2}=0$

其定义在(−1, 1)上，有以下边界条件：

$u\left(-1,t\right)=0,u\left(1,t\right)=0$

令其初始条件为：

$u\left(x,0\right)=\sin \pi x$

已知此方程的精确解为：

$u\left(x,t\right)={\text{e}}^{-{\pi }^{2}t}\sin \pi x$

3.2. 数值求解

利用上述方法先将初始条件展开为有限项切比雪夫多项式的和，初始条件对应时间为0的状态，即：

$\sin \pi x=\underset{k=0}{\overset{N}{{\displaystyle \sum }}}{a}_k\left(0\right)T_k(\; x\; )$

利用Collocation法的对应公式求解得到系数 $a_k\left(0\right)$ ：以x轴上分为32个区间，切比雪夫多项式展开至k = 20，通过matlab程序按照上面的算法首先得到初始条件的切比雪夫系数，u_i为对应sinπx的数值，k = 0时对应的系数计算公式为：

${\hat{u}}_0=\frac{2}{{\bar{c}}_0\ast 32}\underset{i=0}{\overset{32}{{\displaystyle \sum }}}\frac{1}{{\bar{C}}_i}{u}_i\cos 0$

k = 1时对应的值为：

${\hat{u}}_1=\frac{2}{{\bar{c}}_1\ast 32}\underset{i=0}{\overset{32}{{\displaystyle \sum }}}\frac{1}{{\bar{C}}_i}{u}_i\cos \frac{\pi i}{32}$

以此类推可以得到时间为0时所有的系数，再将以k从大到小的顺序先得到一次导数对应的系数，再以相同的方法得到二阶导数的系数，对其求两阶导得到对应的系数 $a_k^{\left(2\right)}\left(0\right)$ ，对应函数 ${u^{\prime }}_t={u^{″}}_x$ ，右侧等于

$\underset{k=0}{\overset{N}{{\displaystyle \sum }}}{a}_k^{\left(2\right)}\left(0\right)T_k(\; x\; )$

此时左侧可以使用数值方法进行估计，若采用向前欧拉法，左侧即为 $\left(a_k\left(1\right)T_k\left(x\right)-a_k\left(0\right)T_k\left(x\right)\right)/\Delta t$ ，令左右两侧相等，可以得到

$a_k\left(1\right)=a_k\left(0\right)+\Delta t{a}_k^{\left(2\right)}(\; 0\; )$

类似的，也可以使用改进欧拉法或龙格库塔法等增加拟合的精度。以时间为依据不断迭代，可以得到a_k关于时间t的变化情况，而通过 $a_k\left(t\right)$ 可以推出函数在不同时刻与x坐标的值 [10] 。

3.3. 计算优化——FFT

使用谱方法计算谱系数会占用大量资源，故在方法出现的若干年间一直没有像有限元法与有限差分法一样得到广泛应用，直到快速傅里叶变换(FFT)的出现，快速傅里叶变换虽然没有为傅里叶变换增加新的理论内容，但其使原先N²复杂度降低到Nlog₂N [11] ，故选取点数时也需要满足为2的整数幂次，快速傅里叶变换使系数计算效率大大提高，谱方法也再次获得广泛应用。

离散傅里叶变换的公式如下：

$\hat{x}\left[k\right]=\underset{n=0}{\overset{N-1}{{\displaystyle \sum }}}{\text{e}}^{-i\frac{2\pi }{N}nk}x\left[n\right]$

将 ${\text{e}}^{-i\text{2}\pi }$ 简记为ω， ${\text{e}}^{-i\frac{2\pi }{N}nk}$ 表示为 ${\omega }_N^{nk}$ ，由其性质我们可以得到1) ${\omega }_N^k={\omega }_N^{k+N}$ 2) ${\omega }_N^k=-{\omega }_N^{k+N/2}$ 3)

${\omega }_N^k={\omega }_{N/2}^{k/2}$ ，若以N = 4为例，原来出现的 ${\omega }_4^i\left(i=1,2,3,\cdots ,9\right)$ 则仅剩下 ${\omega }_4^1$ 与 ${\omega }_2^1$ 两个值，公式的右侧写成矩阵形式并利用这一特点变换有以下：

$\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& {\omega }_4^1& {\omega }_2^1& -{\omega }_4^1\\ 1& {\omega }_2^1& 1& {\omega }_2^1\\ 1& -{\omega }_4^1& {\omega }_2^1& {\omega }_4^1\end{array}\right]\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\\ x_3\end{array}\right]=\left|\begin{array}{c}{x}_0+x_1+x_2+x_3\\ x_0+{\omega }_4^1{x}_1+{\omega }_2^1{x}_2-{\omega }_4^1{x}_3\\ x_0+{\omega }_2^1{x}_1+x_2+{\omega }_2^1{x}_3\\ x_0-{\omega }_4^1{x}_1+{\omega }_2^1{x}_2+{\omega }_4^1{x}_3\end{array}\right|=\left[\begin{array}{c}\left(x_0+x_2\right)+\left(x_1+x_3\right)\\ \left(x_0-x_2\right)+{\omega }_4^1\left(x_1-x_3\right)\\ \left(x_0+x_2\right)-\left(x_1+x_3\right)\\ \left(x_0-x_2\right)-{\omega }_4^1\left(x_1-x_3\right)\end{array}\right]$

可以看出对其奇偶项可以进行拆分并提取公共项

$\{\begin{array}{l}X\left(k\right)=\underset{n=0,1}{{\displaystyle \sum }}{x}_{2n}{\omega }_2^{nk}+{\omega }_4^k\underset{n=0,1}{{\displaystyle \sum }}{x}_{2n+1}{\omega }_2^{nk}k=0,1\\ X\left(k+2\right)=\underset{n=0,1}{{\displaystyle \sum }}{x}_{2n}{\omega }_2^{nk}-{\omega }_4^k\underset{n=0,1}{{\displaystyle \sum }}{x}_{2n+1}{\omega }_2^{nk}k=0,1\end{array}$

在这个过程中，我们仅需计算前2个点，后2个的点利用对称性易得。若N的值为2次幂时，不断拆分奇偶下标项的过程可以一直递归下去，故其计算量为0 (Nlog₂N)。

3.4. 误差分析

3.4.1. 实际误差分析

选取恰当的x步长、展开次数与时间步长使结果随时间收敛，并在运算过程中选取其中部分时间拟合数值与实际数值对比进行分析误差。

T = 0.01

T = 0.05

选取两个时间点的数值解与精确解进行对比，发现整体符合，随着x轴上选取点数的增加，数值解的值也与精确值更加接近。

3.4.2. 理论误差分析

无论是应用Galerkin法还是插值(Collocation)法，近似解与精确解之间都会存在误差，其取决于解的光滑性，一般来说精确解越光滑近似解就越准确。Canuto (1988)与Mercier (1989) [3] 进行过误差理论分析，在 $H_{\omega }^P\left(-1,1\right)$ 范数中，使用Galerkin法误差满足

${\Vert u-u_N\Vert }_{H_{\omega }^P\left(-1,1\right)}\le C{N}^{-\frac{1}{2}+2p-m}{\Vert u\Vert }_{H_{\omega }^P\left(-1,1\right)}$

在同样的范数下，使用Collocation法的误差满足

${\Vert u-u_N\Vert }_{H_{\omega }^P\left(-1,1\right)}\le C{N}^{2p-m}{\Vert u\Vert }_{H_{\omega }^P\left(-1,1\right)}$

同时，在处理N阶切比雪夫近似时，舍入误差(计算值与真实值误差)是必须考虑的问题。一个造成原因是计算切比雪夫展开的导数系数采用了公式计算，其不必借助FFT算法，计算简单，而微分矩阵却可以到达更好的精确度。造成它的一个原因是矩阵的元的维度的巨大差异，如在一阶微分矩阵中，最小的元是 $\text{O}\left(1\right)$ 而最大的是 $\text{O}\left(N^2\right)$ ；另一个原因可以在矩阵元素的计算中找到，它涉及几乎相等的数的减法，Rothman (1991)观察到并提出了一种减小误差的简单方法：使用三角恒等式来表示 $\left(x_i-x_j\right)$ 与 ${\left(1-x_i\right)}^2$ ，令

$x_i-x_j=2\sin \frac{\left(j+i\right)\pi }{2N}\sin \frac{\left(j-i\right)\pi }{2N},i-x_i^2={\sin }^2\frac{i\pi }{N}$

数值实验展现了这一简单过程的有效性 [4] [12] [13] ，同时在二阶微分中使用微分矩阵D²效果甚至更差。另一个舍入误差产生的原因是它不能很准确得展现一个常数的微分计算矩阵 [14] [15] ，即：

$\underset{j=0}{\overset{N}{{\displaystyle \sum }}}{d}_{i,j}^{\left(1\right)}=0,i=0,\cdots ,N$

3.4.3. 总结

整体过程展现谱方法的有效性，提供了偏微分方程区别于有限元法的一种数值解法，方法逻辑顺序清晰理论成熟，类似的方法可以应用于不同的问题，解决各种无法得到数值解的微分方程问题。