在本文中, 将考虑如下形式的Cahn-Hilliard方程^[1]的数值解：

$ \begin{equation} \left \{ \begin{aligned} u_{t}\; \; &= M \Delta W, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; (x, t) \in \Omega \times (0, T], \\ W \; \; & = - \Delta u + \frac{1}{\epsilon^{2}} f(u), \; \; \; \; \; \; \; \; (x, t) \in \Omega \times (0, T], \\ u(x, 0) & = \; \; \; u_{0}, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; x \in \Omega, \\ \frac{\partial u}{\partial n}\; \; &= \frac{\partial W}{\partial n}=0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; x \in \partial \Omega, \end{aligned} \right. \end{equation} $

(1.1)

这里$ \frac{\partial u}{\partial n}=\nabla u \cdot n $, 其中$ n $指的是单位外法向量, $ F(u)=\frac{1}{4} (u^{2}-1)^{2} $是双阱势能, 函数$ f(u)=F'(u)=u^{3}-u $, $ \Omega $是有界区域, $ \epsilon $是界面参数, $ M $为相场迁移率, $ T $表示最终时刻, $ u $指的是相场中某物质的浓度. 其中能量函数为

$ \begin{equation} E(u) = \int_{\Omega} \left ( \frac{1}{2}|\nabla u|^{2} + \frac{1}{\epsilon^{2}}F(u) \right) dx. \end{equation} $

(1.2)

Cahn-Hilliard方程（简记为CH）, 最初是在1958年由Cahn和Hilliard^[1]提出的, 主要是描述二元合金中相分离, 粗化, 合并, 溶解, 迁移等现象, 常应用在材料科学, 相分离等领域. CH方程是一个刚性的, 非线性四阶偏微分方程, 在进行求解时, 涉及到时间步长的选取, 以及需要对高阶空间进行求导, 因此, 往往很难找到解析解. 另外, 由于我们是在一个有界的区域里面考虑CH方程的求解, 所以边界条件的处理对我们的数值求解有着重要的影响, . 如^[7]中提到用有限元法求解基于SAV格式的CH方程的动态边界问题.

数值模拟方法是求解CH方程的一个重要手段之一. 近些年来, 求解CH方程的数值方法有很多, 比如有限差分法^{[2, 3]}, 谱方法^{[4, 5]}, 有限元^{[6, 7]}等. 虽然这些方法在求解CH方程上的研究成果很多, 但仍然面临着许多的挑战. 众所周知, 相分离的过程都满足能量耗散定理, 离散方程得到数值格式时, 需要数值解达到稳定状态, 满足能量耗散定理, 这可能要很大的时间步长^{[8, 9]}. 沈捷等^[8]分别对Allen-Cahn方程和Cahn-Hilliard方程进行数值离散得到数值格式, 并对此数值格式进行修缮, 得到稳定的数值格式. 另外, 基于谱伽辽金方法得到弱形式的基础上分别给出了两种方程的稳定性分析和误差估计. 何银年等^[10]用傅里叶谱方法对空间进行离散, 半隐式方法对时间进行离散；同时, 给出一阶和二阶半隐格式的收敛性分析和稳定性分析, 并通过数值实验验证了较大时间步长方法的有效性. 与经典的Cahn-Hilliard模型不同, 杜强等^[11]提出了非局部Cahn-Hilliard模型, 即具有非局部扩散算子的Cahn-Hilliard方程. 在文献中提到两种不同的数值格式：一阶半隐格式和二阶半隐格式并证明了能量稳定性, 最后给出数值实验验证得到的结果是一致的. 李东等^{[5, 12]}用半隐式傅里叶方法离散二维CH方程得到二阶半隐数值格式, 并提出一种新的稳定方法来证明修缮后的能量函数的稳定性, 并进一步扩展到了三维空间.

谱方法是一种高效的数值方法, 尤其是切比雪夫谱方法^[13-17], 它是求解偏微分方程的一个重要手段之一. 对于很多非周期性问题, 应用傅里叶谱方法^[17]进行求解成本比较高, 而切比雪夫谱方法求解这类问题就十分有优势；通过修改基函数或者导函数以适应各种非周期边界条件. 相比于其他谱方法, 在光滑解的情况下, 切比雪夫谱具有指数收敛性^[13], 即在高精度计算中十分有效, 特别是对应于高阶的偏微分方程, 利用切比雪夫微分矩阵进行相应的计算高阶导数, 这可以避免有限差分中对于高阶导数计算产生的复杂性.

本文主要利用切比雪夫谱方法的优势来求解CH方程. 主要内容如下：第二节, 介绍切比雪夫谱方法的定义, 给出选取不同的类型得到不同的切比雪夫微分矩阵类型. 第三节引出CH方程的半隐格式, 并利用微分矩阵进行离散求解. 第四节, 在第三节的基础上, 增加稳定项, 并利用切比雪夫谱方法对空间进行离散, 时间仍然是半隐式方法离散, 得到数值解. 第五章结合第三、四节, 给出二阶稳定化半隐格式, 与三四节同理, 离散方程得到数值格式. 第六节, 给出数值算例, 比较三种格式的不同及其稳定性情况.

切比雪夫谱方法是谱方法中的一种高效数值算法, 其定义是选取正交多项式为切比雪夫多项式的谱方法即为切比雪夫谱方法.

这里, 我们选取的是第一类切比雪夫多项式, 可表示为

$ \begin{equation} T_{n}(x) = \cos (n \arccos x), \; \; \; \; \; \; n \geq 0, \; \; \; \; \; x \in \Omega, \end{equation} $

(2.1)

通常, 令$ x=\cos \theta $, 则$ \theta = \arccos x $, 那么切比雪夫多项式可简写为$ T_{n}(x)=\cos (n \theta ), n \geq 0, x \in \Omega $. 该多项式满足一些性质, 如：

(1) 递归性

$ \begin{equation} \begin{aligned} T_{n+1}(x)& = 2x T_{n}(x) - T_{n-1}(x) , \; \; \; \; \; \; n \geq 0, \\ T_{0}(x)&=1, \; \; T_{1}(x) =x, \end{aligned} \end{equation} $

(2.2)

与它的导数的递归关系为

$ \begin{equation} 2T_{n}(x) = \frac{1}{n+1} T'_{n+1} (x) - \frac{1}{n-1} T'_{n-1}(x), \; \; \; \; \; \; n \geq 2, \end{equation} $

(2.3)

由余弦函数的和积化差性质可得到

$ \begin{equation} 2T_{m}(x)T_{n}(x) = T_{m+n}(x) - T_{m-n} (x) , \; \; \; \; \; \; n \geq 2. \end{equation} $

(2.4)

(2) 正交性

$ \begin{equation} \begin{aligned} \int_{-1}^{1}T_{n}(x)T_{m}(x)\frac{1}{\sqrt{1-x^{2}}}dx=\frac{c_{n}\pi}{2}\delta_{mn} \quad c_{0}=2, \quad c_{n}=1(n\geq1), \end{aligned} \end{equation} $

(2.5)

其中

$ \begin{equation} \delta_{mn} =\left\{ \begin{aligned} &0, \; \; \; \; \; \; \; \; \; \quad m\neq n, \\ &1, \; \; \; \; \; \; \; \; \; \quad m=n, \end{aligned} \right. \end{equation} $

(2.6)

权重函数为

$ \begin{equation} w=\frac{1}{\sqrt{1-x^{2}}}. \end{equation} $

(2.7)

而对于离散点的选取有三种选择：

(1) Chebyshev-Guass (CG) 点

$ \begin{equation*} x_{j}=\cos(\frac{(2j+1)\pi}{2N+2}), \quad w_{j}=\frac{\pi}{N+1}, \quad0\leq j\leq N. \end{equation*} $

(2) Chebyshev-Guass-Radau (GCR) 点

$ \begin{equation*} x_{0}=1, w_{0}=\frac{\pi}{2N+1}, x_{j}=\cos(\frac{(2\pi j+1)\pi}{2N+1}), \quad w_{j}=\frac{2\pi}{2N+1}, \quad1\leq j\leq N. \end{equation*} $

(3) Chebyshev-Guass-Lobatto (CGL) 点

$ \begin{equation*} x_{0}=1, x_{N}=-1, w_{0}=w_{N}=\frac{\pi}{2N}, x_{j}=\cos(\frac{(\pi j)\pi}{N}), w_{j}=\frac{\pi}{N}, \quad1\leq j\leq N-1. \end{equation*} $

一般情况下, 对于不同的边界条件问题, 离散点的选取通常不同. 例如求解椭圆和抛物类问题, 常选取CGL点, 然而对于一阶双曲问题常选用CG点和CGR点.

在用切比雪夫谱方法求解时, 常需要构造切比雪夫微分矩阵. 例如, 选取CG型离散点, 并设离散点为$ \{x\}_{j=0}^{N} $. 不妨设$ u\in P_{N} $, 其中$ P_{N} $为$ N $次线性多项式, $ \{L_j(x)\}_{j=0}^N $为拉格朗日基多项式, 则

$ \begin{equation} u^{(m)}=D^{m}u, \quad m\geq1, \end{equation} $

(2.8)

其中$ D $即为切比雪夫微分矩阵, 表达式如下：

$ \begin{equation} \begin{aligned} D & = (d_{kj} = L'_{j}(x_{k}))_{0 \leq k, j \leq N}, \\ u^{(m)} & = (u^{(m)} (x_{0}), \dots, u^{(m)} (x_{N})), \\ u &= u^{0}. \end{aligned} \end{equation} $

(2.9)

对于选取不同的离散点, 微分矩阵的构造不同, 这里只考虑CGL$ (x_{0}=-1, x_{N}=1) $型. 下面是微分矩阵的每个元素的表达式：

$ \begin{equation} \left. d_{kj}=\left\{ \begin{array} {cc}{-\frac{2N^{2}+1}{6}, } & {k=j=0}, \\ {\frac{\widetilde{c}_{k}}{\widetilde{c}_{j}}\frac{(-1)^{k+j}}{x_{k}-x_{j}}, } & \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; {k\neq j, 0\leq k, j\leq N}, \\ {-\frac{x_{k}}{2(1-x_{k}^{2})}, } & \; \; \; \; \; \; \; \; \; \; \; \; \; \; {1\leq k=j\leq N-1}, \\ {\frac{2N^{2}+1}{6}, } & {k=j=N.} \end{array}\right.\right. \end{equation} $

(2.10)

注：$ \widetilde{c}_{0}=\widetilde{c}_{N}=2 $, $ \widetilde{c}_{j}=1, \; 1 \leq j \leq N-1. $

半隐格式是数值求解方法中常用的一种, 在偏微分方程对时间进行离散时, 对等式右边线性部分采用隐式处理, 非线性部分采用显示处理, 这样做的优势相比完全显性格式来说可以显著地减少计算量, 节约成本. 本节考虑采用半隐格式求解具有狄利克雷边界条件和齐次冯诺依曼边界条件的Cahn-Hilliard方程.

$ \begin{equation} \left\{ \begin{array} {l}\frac{u^{n+1}-u^n}{\delta t}=M\Delta\left[-\Delta u^{n+1}+\frac{f(u^{n})}{\epsilon^2}\right], \\ u(x, 0)=u_0, \end{array} \right. \end{equation} $

(3.1)

设其初值条件为

$ \begin{equation} u(x, 0)=cos^{2}(\pi x), \end{equation} $

(3.2)

以及狄利克雷边界, 齐次冯诺依曼边界条件

$ \begin{equation} \begin{aligned} u(-1)&=u(1)=0, \\ u^{\prime}(-1)&=u^{\prime}(1)=0, \end{aligned} \end{equation} $

(3.3)

其中$ {{u}^{n+1}} $表示在时间$ t=t_{n+1} $时刻$ u(x, {{t}_{n+1}}) $的值, $ t={{t}_{n+1}} $是将时间剖分为$ N $等份取$ n+1 $时刻的值：$ 0=t_{0}<t_{1}<t_{2}<\cdots<t_{N} $, 其中$ \delta t:=\frac{{{t}_{N}}}{N} $表示时间步长, $ u_{0} $表示在时间为0时刻的初始值.

首先, 我们用切比雪夫谱方法离散空间变量：

$ \begin{equation} \frac{u^{n+1}-u^{n}}{\delta t}=-MD^{4}u^{n+1}+MD^{2}\frac{f(u^{n})}{\epsilon^{2}}, \end{equation} $

(3.4)

整理得

$ \begin{equation} u^{n+1}-u^{n}=-M\delta tD^{4}u^{n+1}+M\delta tD^{2}\frac{f(u^{n})}{\epsilon^{2}}, \end{equation} $

(3.5)

合并同类项得

$ \begin{equation} (I+M\delta tD^{4})u^{n+1}=(I-\frac{M\delta tD^{2}}{\epsilon^{2}})u^{n}+\frac{M\delta tD^{2}}{\epsilon^{2}}(u^{n})^{3}. \end{equation} $

(3.6)

为简便运算, 记

$ \begin{equation} \left\{ \begin{aligned} A &=I+M\delta tD^{4}, \\ B_{1}&=I-\frac{M\delta tD^{2}}{\epsilon^{2}}, \\ B_{2}&=\frac{M\delta tD^{2}}{\epsilon^{2}}, \\ b &= B_{1}u^{n} +B_{2} (u^{n})^{3}, \end{aligned} \right. \end{equation} $

(3.7)

其中, $ I $表示与空间相同阶的单位矩阵, $ D $是切比雪夫微分矩阵. 则我们的离散方程最终就化简为

$ \begin{equation} Au^{n+1}=b. \end{equation} $

(3.8)

那么, 只需求解上述线性方程组, 即可得到半隐格式的数值解.

本节考虑Cahn-Hilliard方程的一阶稳定化半隐数值格式. 此格式在是第三节半稳定格式的右边添加稳定项$ \frac{ \Delta S(u^{n+1}-u^n)}{\epsilon^2} $, 其中$ S $是稳定化常数, 从而得到一阶稳定化半隐格式^[18]（记为一阶SSI1）.

首先, 在一维区间$ [-1, 1] $上考虑如下格式的Cahn-Hilliard方程：

$ \begin{equation} \left\{ \begin{array} {ll} \frac{u^{n+1}-u^n}{\delta t}&=M\Delta W^{n+1}, \\ W^{n+1} & =-\Delta u^{n+1}+\frac{S}{\epsilon^2}(u^{n+1}-u^n)+\frac{1}{\epsilon^2}f(u^n), \\ u(x, 0)&=cos^{2}(\pi x), \\ u^{\prime}(-1)&=u^{\prime}(1)=0, \\ u(-1)&=u(1)=0, \end{array}\right. \end{equation} $

(4.1)

其中, $ f(u^n)=(u^n)^3-u^n $. 上述等式等价为：

$ \begin{equation} \frac{u^{n+1}-u^{n}}{\delta t}=M\Delta[-\Delta u^{n+1}+\frac{S}{\epsilon^{2}}(u^{n+1}-u^{n})+\frac{1}{\epsilon^{2}}f(u^{n})]. \end{equation} $

(4.2)

首先对空间离散得：

$ \begin{equation} \frac{u^{n+1}-u^{n}}{\delta t}=MD^{2}[-D^{2}u^{n+1}+\frac{S}{\epsilon^{2}}(u^{n+1}-u^{n})+\frac{1}{\epsilon^{2}}f(u^{n})], \end{equation} $

(4.3)

进一步化简得

$ \begin{equation} [1+MD^{4}\delta t-\frac{MD^{2}S\delta t}{\epsilon^{2}}]u^{n+1}=[I-\frac{MD^{2}S\delta t}{\epsilon^{2}}]+\frac{\delta tMD^{2}}{\epsilon^{2}} f(u^{n}), \end{equation} $

(4.4)

并将$ f(u^{n}) $带入, 最终化简得

$ \begin{equation} [I+MD^{4}\delta t-\frac{MD^{2}S\delta t}{\epsilon^{2}}]u^{n+1}=[I-\frac{MD^{2}S\delta t}{\epsilon^{2}}-\frac{\delta tMD^{2}}{\epsilon^{2}}]u^{n}+\frac{\delta tMD^{2}}{\epsilon^{2}}(u^{n})^{3}, \end{equation} $

(4.5)

其中$ I $为阶单位矩阵, 为了后续计算方便, 记

$ \begin{equation} \left \{ \begin{aligned} & U=I-\frac{MD^{2}S\delta t}{\epsilon^{2}}, \\ & A=U+M\delta tD^{4}, \\ & b=(U-\frac{M\delta tD^{2}}{\epsilon^{2}})u^{n}+\frac{M\delta tD^{2}}{\epsilon^{2}}(u^{n})^{3}, \end{aligned} \right. \end{equation} $

(4.6)

则上述的一阶格式就化简为线性方程组

$ \begin{equation} Au^{n+1}=b. \end{equation} $

(4.7)

在文献[5]中提到, 当$ S $满足以下条件

$ \begin{equation} S \geq \max\limits_{x \in \Omega} \{ \frac{1}{2} |u^{n}(x)|^{2} + \frac{1}{4} |u^{n+1}(x) + u^{n}(x)|^{2} \} -\frac{1}{2}, \end{equation} $

(4.9)

能量是稳定的.

与一阶稳定化半隐格式的构造思想是类似的, 在二阶线性半隐格式（记为SSI2）^[13]上添加稳定项, 且它的截断误差阶为, 从而得到二阶SSI格式. 此格式是由^[10]在求解二维的MBE方程时提出的.

本节采用二阶向后微分公式(BDF2)^[19]逼近该数值格式. 在1952年, BDF方法最初由C. F. Curtiss和J. O. Hirschfelder提出^[20], 用于求解具有刚性的微分方程, 即解决刚性问题. 并在1971年, 由Gear^[21]在其经典著作中研究了BDF(包括BDF1到BDF6的所有公式及性质)的数学分析和稳定性, 这个方法是如今数值求解刚性问题的标准之一.

现在来看BDF的公式, 考虑ODE方程

$ \begin{equation} \frac{dy}{dt} = f(t, y), \end{equation} $

(5.1)

则运用BDF公式进行隐式离散后

$ \begin{equation} \frac{y^{n+1} - y^{n}}{\delta t}= f(t_{n}, y^{n+1}), \end{equation} $

(5.2)

那么在此基础上运用BDF2公式, 则离散形式为：

$ \begin{equation} \frac{3y^{n+2} - 4y^{n+1} + y^{n}}{2 \delta t} = f(t_{n+2}, y_{n+2}), \end{equation} $

(5.3)

其中, $ \delta t $表示时间步长, 计算BDF2时, $ y^{0}, y^{1} $需要通过BDF公式计算得来, 再带入BDF2公式进行迭代运算.

采用BDF2公式求解的Cahn-Hilliard方程的二阶SSI格式^[15]为

$ \begin{equation} \left\{ \begin{aligned} &\frac{3u^{n+1}-4u^n+u^{n-1}}{2\delta t} =M\Delta W^{n+1} , \\ &W^{n+1} =-\Delta u^{n+1}+\frac{2f(u^n)-f(u^{n-1})}{\epsilon^2} +\frac{S(u^{n+1}-2u^n+u^{n-1})}{\epsilon^2}, \end{aligned} \right. \end{equation} $

(5.4)

同样设初始值为

$ \begin{equation} u(x, 0)=cos^{2}(\pi x), \end{equation} $

(5.5)

和Dirichlet, Neumann条件

$ \begin{equation} \begin{aligned} u(-1) & =u(1)=0, \\ u^{\prime}(-1) & =u^{\prime}(1)=0, \end{aligned} \end{equation} $

(5.6)

其中与一阶SSI一样, $ f(u)=u^{3}-u $, $ S $是稳定化常数, $ \delta t $是时间步长. 上述二阶SSI格式(5.4)式等价为

$ \begin{equation} \begin{aligned} \frac{3u^{n+1}-4u^{n}+u^{n-1}}{2\delta t} =M\Delta[-\Delta u^{n+1}+\frac{2f(u^{n})-f(u^{n-1})}{\epsilon^{2}} +\frac{S(u^{n+1}-2u^{n}+u^{n-1})}{\epsilon^{2}}]. \end{aligned} \end{equation} $

(5.7)

现用切比雪夫谱方法求解该二阶格式(5.4), 对空间进行离散得到

$ \begin{equation} \begin{aligned} \frac{3u^{n+1}-4u^{n}+u^{n-1}}{2\delta t}=MD^{2}[-D^{2}u^{n+1}+\frac{2f(u^{n})-f(u^{n-1})}{\epsilon^{2}} +\frac{S(u^{n+1}-2u^{n}+u^{n-1})}{\epsilon^{2}}], \end{aligned} \end{equation} $

(5.8)

进一步化简得

$ \begin{equation} \begin{aligned} \begin{aligned} 3u^{n+1}-4u^{n}+u^{n-1} \end{aligned} =&-2\delta tMD^{4}u^{n+1}+\frac{2\delta tMD^{2}S}{\epsilon^{2}} u^{n+1}\\ & +2\delta tMD^2\frac{2(u^n)^3-(u^{n-1})^3-2(1+S)u^n+(1+S)u^{n-1}}{\epsilon^2}, \end{aligned} \end{equation} $

(5.9)

最后得到

$ \begin{equation} \begin{aligned} &(3I+2\delta tMD^{4}-\frac{2\delta tMSD^{2}}{\epsilon^{2}})u^{n+1}\\ =&(4I-\frac{4\delta tM(1+S)D^2}{\epsilon^{2}})u^n +(\frac{2\delta tM(1+S)D^{2}}{\epsilon^{2}}-I)u^{n-1} +2\delta tMD^2\frac{2(u^n)^3-(u^{n-1})^3}{\epsilon^2}. \end{aligned} \end{equation} $

(5.10)

为了简便运算, 记

$ \begin{equation} \left \{ \begin{aligned} & U_{1}=\frac{2\delta tMD^{2}}{\epsilon^{2}}, \\ & U_{2}=I-\frac{2\delta tM(1+S)D^{2}}{\epsilon^{2}}, \\ & A=3I+2\delta tMD^{4}-U_{1}, \\ & B_{1}=2I+2U_{2}, \\ & B_{2}=-U_{2}, \\ & B_{3}=\frac{U_{1}}{S}, \\ & b=B_{1}u^{n}+B_{2}u^{n-1}+B_{3}(2(u^{n})^{3}-(u^{n-1})^{3}), \end{aligned} \right. \end{equation} $

(5.11)

因此, 可将二阶SSI格式记为如下

$ \begin{equation} A u^{n+1}=b. \end{equation} $

(5.12)

在进行迭代运算求解时, $ u^{0}, u^{1} $的值需要通过一阶SSI格式计算得出, 再带入二阶格式进行运算, 直到达到最后运算时间时, 即可得出近似解.

为方便接下来的计算过程, 在此引入一些符号：$ \delta_{t} u^{n+1} := u^{n+1}-u^{n} $, $ \delta_{tt} u^{n+1}:=u^{n+1}-2u^{n} + u^{n-1} $, $ D_{\delta t} u^{n+1}:=\frac{3u^{n+1} - 4u^{n}- u^{n-1}}{2 \delta t} $. 对于前一节提到的二阶稳定化数值格式可化为如下(5.13)式：

$ \begin{equation} \left \{ \begin{aligned} D_{\delta t} u^{n+1} & = M \Delta W^{n+1}, \\ W^{n+1} & = - \Delta u^{n+1} + \frac{1}{\epsilon^{2}} (2f(u^{n})-f(u^{n-1})) + \frac{S}{\epsilon^{2}} \delta_{tt} u^{n+1}. \end{aligned} \right. \end{equation} $

(5.13)

对于任意$ u \in L^{2}_{0} (\Omega) $, 令$ -\Delta^{-1} u : = u_{1} \in H^{1}(\Omega) \bigcap L^{2}_{0} (\Omega) $, 其中$ u_{1} $是下面具有冯诺依曼边界条件的方程的解：

$ \begin{equation*} \left \{ \begin{aligned} - \Delta u_{1} & = u , \\ \frac{\partial u_{1}}{\partial n} & =0, \end{aligned} \right. \end{equation*} $

其中定义范数$ \| u\|_{-1} := \sqrt {(u, -\Delta^{-1} u)} $.

为证明数值格式的能量稳定性, 以及在后续证明中对于$ f(u) $进行Taylor展开并化简时需要对它的导数进行处理, 所以假设方程(1.1)中的函数$ f'(u) $是要满足一致有界的^[22], 即：

$ \begin{equation} \max\limits_{u \in \mathbf{R}} |f'(u)| \leq L, \end{equation} $

(5.14)

其中$ L $是非负常数. 并且在文献[21]中证明了二阶稳定化半隐格式是满足

$ \begin{equation} \frac{1}{|\Omega|} \int u^{0} dx = \frac{1}{|\Omega|} \int u^{1} dx = m_{0}, \end{equation} $

(5.15)

从而由数学归纳法得到

$ \begin{equation} \frac{1}{|\Omega|} \int u^{n+1} dx = m_{0}, \; \; \; n=1, 2, \dots, N, \end{equation} $

(5.16)

则$ \delta_{t} u^{n+1} \in L^{2}_{0}, \; \; n=0, \dots, N. $

定理5.1  在满足(5.15)式以及(5.16)式的情况下, %当$ S<L $时, 二阶SSI数值格式(5.13) 满足如下的能量稳定性

$ \begin{equation} \begin{aligned} E_{S}^{n+1} & \leq E_{S}^{n} - \frac{1}{4M\delta t } \|\delta_{tt}u^{n+1}\|^{2}_{-1} - \frac{1}{2} \|\ \nabla \delta_{t} u^{n+1}\|^{2}\\ &-\frac{1}{M\delta t} \| \delta_{t} u^{n+1}\|^{2}_{-1} + \frac{L}{\epsilon^{2}} \| \delta_{t} u^{n+1}\|^{2} - \frac{S}{2\epsilon^{2}} \| \delta_{tt} u^{n+1}\|^{2}, \end{aligned} \end{equation} $

(5.17)

其中

$ \begin{equation} E_{S}^{n+1} = E(u^{n+1}) + \frac{1}{4M\delta t } \|\delta_{t}u^{n+1}\|^{2}_{-1} + (\frac{L}{2\epsilon^{2}} + \frac{S}{2\epsilon^{2}}) \| \delta_{t} u^{n+1}\|^{2}. \end{equation} $

(5.18)

证  在(5.13)式左右两边分别与$ \frac{(-\Delta)^{-1} \delta_{t} u^{n+1}}{M} $和$ -\delta u^{n+1} $做内积, 再利用分布积分公式, 散度定理, 和以下关系式：

$ \begin{equation*} \begin{aligned} (a, a-b) & = \frac{1}{2} (\|a\|^{2} - \|b\|^{2} + \|a - b\|^{2}), \\ (a, b)\; \; \; \; & = \frac{1}{2} (\|a\|^{2} + \|b\|^{2} - \|a - b\|^{2}), \end{aligned} \end{equation*} $

得到

$ \begin{equation} \begin{aligned} (\frac{1}{M}D_{\delta t} u^{n+1}, (-\Delta)^{-1} \delta u^{n+1}) &= -(\nabla u^{n+1}, \nabla \delta_{t} u^{n+1}) - \frac{1}{\epsilon^{2}} (2f(u^{n})-f(u^{n-1}), \delta_{t} u^{n+1}) \\ &-\frac{S}{\epsilon^{2}} (\delta_{tt} u^{n+1}, \delta_{t} u^{n+1}). \end{aligned} \end{equation} $

(5.19)

再做简单的变换得到

$ \begin{equation} I_{1}+I_{2}+I_{3}+I_{4}=0, \end{equation} $

(5.20)

其中

$ \begin{equation} \begin{aligned} I_{1}&= -(\frac{1}{M}D_{\delta t} u^{n+1}, (-\Delta)^{-1} \delta u^{n+1}) \\ &= -(\frac{1}{M \delta t} \delta_{t} u^{n+1} + \frac{1}{2 M \delta t} \delta_{tt} u^{n+1}, (-\Delta)^{-1} \delta_{t} u^{n+1})\\ &= -\frac{1}{M \delta t} \| \delta_{t} u^{n+1}\|_{-1}^{2} -\frac{1}{4 M \delta t} (\| \delta_{t} u^{n+1}\|_{-1}^{2} - \| \delta_{t} u^{n}\|_{-1}^{2} + \| \delta_{tt} u^{n+1}\|_{-1}^{2}). \end{aligned} \end{equation} $

(5.21)

$ I_{2} $的表达式如下(5.22)式：

$ \begin{equation} \begin{aligned} I_{2}&= -(\nabla u^{n+1}, \nabla \delta_{t} u^{n+1}) \\ &= -\frac{1}{2} (\| \nabla u^{n+1}\|^{2} - \| \nabla u^{n}\|^{2} + \| \nabla \delta_{t} u^{n+1}\|^{2}). \end{aligned} \end{equation} $

(5.22)

$ I_{3} $为：

$ \begin{equation} \begin{aligned} I_{3}&= -\frac{S}{\epsilon^{2}}(\delta_{tt} u^{n+1}, \delta_{t} u^{n+1}) \\ &= -\frac{S}{2 \epsilon^{2}} \| \delta_{t} u^{n+1}\|^{2} \| + \frac{S}{2 \epsilon^{2}} \| \delta_{t} u^{n}\|^{2} \| - \frac{S}{2 \epsilon^{2}} + \| \delta_{tt} u^{n+1}\|^{2}. \end{aligned} \end{equation} $

(5.23)

求解$ I_{4} = (2f(u^{n}) - f(u^{n-1}), \delta_{t} u^{n+1}) $, 利用Taylor展开公式在$ u^{n} $处展开得到

$ \begin{equation} F(u^{n+1}) = F(u^{n}) + f(u^{n})(u^{n+1}-u^{n}) + \frac{f'(\xi^{n})}{2} (u^{n+1}-u^{n})^{2}, \end{equation} $

(5.24)

则可得到

$ \begin{equation} (u^{n+1}-u^{n}, f(u^{n}))=(1, F(u^{n+1})-F(u^{n})) - \frac{f'(\xi^{n})}{2} (u^{n+1}-u^{n})^{2}, \end{equation} $

(5.25)

从而根据上(5.25)式可得到

$ \begin{equation} (u^{n+1}-u^{n}, f(u^{n})-f(u^{n-1}))=(u^{n+1}-u^{n}, f'(\xi^{n-1})(u^{n}-u^{n-1})), \end{equation} $

(5.26)

则$ I_{4} $为

$ \begin{equation} \begin{aligned} I_{4} &= (f(u^{n}), \delta_{t} u^{n+1}) + (f(u^{n}) - f(u^{n-1}), \delta_{t} u^{n+1}) \\ &= (f(u^{n}), u^{n+1}-u^{n}) + (f(u^{n}) - f(u^{n-1}), u^{n+1}-u^{n}). \end{aligned} \end{equation} $

(5.27)

其中$ \xi^{n} $是介于$ u^{n} $与$ u^{n+1} $之间, $ \xi^{n-1} $是介于$ u^{n-1} $与$ u^{n} $之间. 再对(5.26), (5.27)两边分别同时乘以$ \frac{1}{\epsilon^{2}} $, 并结合(5.20)-(5.27) 式得到

$ \begin{equation} \begin{aligned} &\frac{1}{\epsilon^{2}}(F(u^{n+1})-F(u^{n}), 1) + \frac{1}{2} (\| \nabla u^{n+1} \|^{2}-\| \nabla u^{n}\|^{2}) \\ & + \frac{1}{4 M\delta t} (\| \delta t u^{n+1}\|^{2}_{-1} - \| \delta t u^{n}\|^{2}_{-1}) + (\frac{L}{2\epsilon^{2}} + \frac{S}{2\epsilon^{2}}) (\|\delta_{t} u^{n+1}\|^{2}- \|\delta_{t} u^{n}\|^{2})\\ \leq& -\frac{1}{4M\delta t} \| \delta_{tt} u^{n+1}\|^{2}_{-1} - \frac{1}{M\delta t} \| \delta_{t} u^{n+1}\|^{2}_{-1} - \frac{1}{2} \| \nabla \delta_{t} u^{n+1}\|^{2}\\ &+ \frac{3L}{2\epsilon^{2}} \| \delta_{t} u^{n+1}\|^{2} %- \frac{S}{2\epsilon^{2}} \| \delta_{tt} u^{n+1}\|^{2} + (\frac{L}{2\epsilon^{2}} - \frac{S}{2\epsilon^{2}}) \| \delta_{tt} u^{n+1}\|^{2}. \end{aligned} \end{equation} $

(5.28)

证毕.

在本节, 将要通过数值算例取求解三种数值格式的数值解的稳定性. 在下面的例子中考虑的空间都是一维空间$ \Omega=[-1, 1] $, $ N $是空间离散网格数, $ x_{i}=cos(\frac{(i-1)\theta}{N-1}) $是网格点. 对时间区间$ [0, T] $进行划分：$ 0=t_{0}<t_{1}<t_{2}<\cdots<t_{N_{1}}=T $, 其中时间步长$ \delta t=T/N_{1} $, 根据(4.8) 式和定理5.1, 我们可选取稳定化常数$ S $为$ S=0, 1, 2 $, $ S=0 $即为半隐格式, 也即非稳定格式.

(1). 半隐格式

设置各参数分别为, $ N=128, M=0.01, \epsilon=0.1 $, $ \delta t=10^{-4} $. 当时间到达$ T=2 $时, 数值解的稳定性情况如图 1所示.

图 1中的(a), (b)图显示的初始值是$ T=0 $时的, 瞬态解是$ T=0.06 $时的, 稳态解是$ T=0.2 $时的, (c)图显示的是$ T=0.2 $时的等高线图. 从图 1中可以看到, 当$ T=2 $时, (3.1)式的数值解达到稳定, CPU耗时47.69（注：参数采用的是无量纲化参数）.

(2). 一阶稳定化格式

分别考虑参数$ N=150 $, $ M=0.01 $, $ \epsilon=0.02 $, $ S=2 $, $ \delta = 2 \times 10^{-5} $. 采用一阶稳定化半隐递推格式直至时间到$ T=0.2 $, 数值解的稳定情况如图 2所示：

从图 2中能看到, 当时间到达$ T=0.2 $时, (4.1)式的数值解达到稳定状态.

现在要研究稳定性参数S给数值解精度带来的影响, 令$ N=200 $, 并固定参数$ M=0.01 $, $ \epsilon=0.1 $, 最后, 画出$ T=0.04 $和$ T=5 $时, 对应不同时间步长的图.

观察图 3可得到：对于左边的瞬态解, 即当取时间$ T=0.04 $时, 取时间步长一致为$ \delta t =10^{-4} $, 稳定化$ S=0, 1, 2 $参数有相似的数值解, 而当时间取$ T=5 $时, 三种稳定化参数对应的数值解高度稳合, 但是时间步长却存在相当大的差异, 稳定化参数$ (S=1, 2) $对应的稳定格式是非稳定格式$ (S=0) $的1000倍.

为了研究数值解的稳定性, 我们选取数值解稳定时的最大时间步长$ \delta t_{max} $进行判断, 即此时若时间步长$ \delta t > \delta t_{max} $时, 数值解将会出现爆破的现象. 例如, 考虑当$ n=256 $, $ \epsilon=0.1 $, $ M = 0.1, 0.01 $时, 对于不同的稳定化参数$ S=0, 1, 2 $, 通过数次实验后, 对应的最大时间如表 1所示：

从表 1中可知, 当数值解达到稳定时, 稳定化格式$ S=1, 2 $的最大时间步长远比非稳定化格式$ S=0 $要大.

(3). 二阶稳定化格式

针对二阶稳定化数值格式, 我们取$ N=150 $, $ M=0.01 $, $ \epsilon=0.02 $, $ S=2 $, $ \delta t =2 \times 10^{-5} $, 采用(37)的格式, 计算直到时间为$ T=0.2 $时, 得到的图如图 4所示.

根据图 4可知, 当时间趋近于$ T=0.2 $时, 数值解逐渐变得稳定.

与一阶格式同理, 参数设置相同, 所得图形如图 5所示：

观察图 5可知, 在时间$ t=0.04 $时所得的瞬时解图中, 当三种不同稳定化参数对应的稳定化数值解十分相似, 几乎不能用肉眼做区分时, 发现时间步长的基本一致. 在时间$ T=5 $时所得的稳态解图中, 当三种不同稳定化常数对应的数值解高度重合时, 发现稳定格式$ (S=1, 2) $的时间步长是非稳定格式$ (S=0) $的1000倍.

稳定性：与一阶格式类似, 选取$ n=256 $, $ \epsilon=0.1 $, $ M = 0.1, 0.01 $时, 对于不同的稳定化参数$ S=0, 1, 2 $, 二阶稳定化格式对应的最大时间如表 2所示.

从表 2中可观察到, 当稳定化参数$ S=1, 2 $数值解稳定时, 最大时间步长大于1, 而对于表中选取的参数, 非稳定性格式($ S=0 $)找不到最大时间步长使得解稳定.