偏微分方程算法之二维初边值问题（交替方向隐（ADI）格式）

一、研究对象

以二维抛物型方程初边值问题为研究对象：

$\left\{\begin{matrix} \frac{\partial u(x,y,t)}{\partial t}-(\frac{\partial^{2}u(x,y,t)}{\partial x^{2}}+\frac{\partial^{2}u(x,y,t)}{\partial y^{2}})=f(x,y,t),(x,y)\in \Omega=[0,a]\times [0,b],0<t\leqslant T,\\ u(x,y,0)=\varphi(x,y),(x,y)\in \Omega,\space\space\space\space(1)\\ u(0,y,t)=g_{1}(y,t),u(a,y,t)=g_{2}(y,t),0\leqslant y\leqslant b,0<t\leqslant T,\\ u(x,0,t)=g_{3}(x,t),u(x,b,t)=g_{4}(x,t),0\leqslant x\leqslant a,0<t\leqslant T \end{matrix}\right.$

为了确保连续性，公式（1）中的相关函数满足：

$g_{1}(0,t)=g_{3}(0,t),g_1(b,t)=g_4(0,t),$

$g_{2}(0,t)=g_{3}(a,t),g_{2}(b,t)=g_{4}(a,t)$

$\varphi(0,y)=g_{1}(y,0),\varphi(a,y)=g_{2}(y,0),$

$\varphi(x,0)=g_{3}(x,0),\varphi(x,b)=g_{4}(x,0)$

二、理论推导

2.1 向前欧拉格式

首先进行网格剖分。将三维长方体空间（二维位置平面+一维时间轴）进行剖分，将区域[0,a]等分m份，区域[0,b]等分n份，区域[0,T]等分 $l$ 份，有：

$x_{i}=i\cdot \Delta x=\frac{ia}{m},0\leqslant i\leqslant m,$

$y_{j}=j\cdot \Delta y=\frac{jb}{n},0\leqslant j\leqslant n,t_{k}=k\cdot \Delta t=\frac{kT}{l},0\leqslant k\leqslant l$

得到网格节点坐标 $(x_{i},y_{j},t_{k})$ 。利用数值法求的数值解 $u(x,y,t)$ 在网格节点 $(x_{i},y_{j},t_{k})$ 处的近似值 $u^{k}_{i,j}$ 。

然后将原方程弱化为节点离散方程，即

$\left\{\begin{matrix} \frac{\partial u}{\partial t}|_{(x_{i},y_{j},t_{k})}- (\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j},t_{k})}=f(x_{i},y_{j},t_{k}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,0<k\leqslant l,\\ u(x_{i},y_{j},t_{0})=\varphi(x_{i},y_{j}),0\leqslant i \leqslant m,0\leqslant j\leqslant n, \space\space\space\space (2)\\ u(x_{0},y_{j},t_{k})=g_{1}(y_{j},t_{k}),u(x_{m},y_{j},t_{k})=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u(x_{i},y_{0},t_{k})=g_{3}(x_{i},t_{k}),u(x_{i},y_{n},t_{k})=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.$

然后利用差商代替微商，取

$\frac{\partial u}{\partial t}|_{(x_{i},y_{j},t_{k})}\approx \frac{u(x_{i},y_{j},t_{k+1})-u(x_{i},y_{j},t_{k})}{\Delta t}$

$\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k})}\approx \frac{u(x_{i+1},y_{j},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i-1},y_{j},t_{k})}{\Delta x^{2}}$ 及

$\frac{\partial^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k})}\approx \frac{u(x_{i},y_{j+1},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j-1},t_{k})}{\Delta y^{2}}$

将上面三式代入公式（2），用数值解代替精确解并忽略高阶项，可得数值格式为

$\left\{\begin{matrix} \frac{u^{k+1}_{i,j}-u^{k}_{i,j}}{\Delta t}-(\frac{u^{k}_{i+1,j}-2u^{k}_{i,j}+u^{k}_{i-1,j}}{\Delta x^{2}}+\frac{u^{k}_{i,j+1}-2u^{k}_{i,j}+u^{k}_{i,j-1}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,0\leqslant k\leqslant l-1,\\ u^{0}_{i,j}=\varphi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\\ u^{k}_{0,j}=g_{1}(y_{j},t_{k}),u^{k}_{m,j}=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u^{k}_{i,0}=g_{3}(x_{i},t_{k}),u^{k}_{i,n}=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.$

记 $r_{1}=\frac{\Delta t}{\Delta x^{2}},r_{2}=\frac{\Delta t}{\Delta y^{2}}$ ，则上式可整理为

$\left\{\begin{matrix} u^{k+1}_{i,j}=r_{1}u^{k}_{i-1,j}+r_{1}u^{k}_{i+1,j}+r_{2}u^{k}_{i,j-1}+r_{2}u^{k}_{i,j+1}+(1-2r_{1}-2r_{2})u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t,1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,0\leqslant k\leqslant l-1,\\ u^{0}_{i,j}=\varphi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\space\space\space\space(3)\\ u^{k}_{0,j}=g_{1}(y_{j},t_{k}),u^{k}_{m,j}=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l, \\ u^{k}_{i,0}=g_{3}(x_{i},t_{k}),u^{k}_{i,n}=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.$

公式（3）的截断误差为 $O(\Delta t +\Delta x^{2}+\Delta y^{2})$ ，可证明其稳定性条件为 $r_{1}+r_{2}\leqslant \frac{1}{2}$ 。可见，时间、空间步长的选择条件苛刻，应用价值不大，需要通过其它方式来降低使用局限性。

2.2 Crank-Nicolson格式

显格式方法稳定性条件差，尝试采用隐格式来计算。将原节点离散方程公式（2）改为：

$\left\{\begin{matrix} \frac{\partial u}{\partial t}|_{(x_{i},y_{j},t_{k+\frac{1}{2}})}- (\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j},t_{k+\frac{1}{2}})}=f(x_{i},y_{j},t_{k+\frac{1}{2}}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,0<k\leqslant l,\\ u(x_{i},y_{j},t_{0})=\varphi(x_{i},y_{j}),0\leqslant i \leqslant m,0\leqslant j\leqslant n, \space\space\space\space (4)\\ u(x_{0},y_{j},t_{k})=g_{1}(y_{j},t_{k}),u(x_{m},y_{j},t_{k})=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u(x_{i},y_{0},t_{k})=g_{3}(x_{i},t_{k}),u(x_{i},y_{n},t_{k})=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.$

然后利用差商代替微商，取

$\frac{\partial u}{\partial t}|_{(x_{i},y_{j},t_{k+\frac{1}{2}})}\approx \frac{u(x_{i},y_{j},t_{k+1})-u(x_{i},y_{j},t_{k})}{\Delta t}$

$\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k+\frac{1}{2}})}\approx \frac{1}{2}[\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k})}+\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k+1})}]\approx\frac{1}{2}[\frac{u(x_{i+1},y_{j},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i-1},y_{j},t_{k})}{\Delta x^{2}}+\frac{u(x_{i+1},y_{j},t_{k+1})-2u(x_{i},y_{j},t_{k+1})+u(x_{i-1},y_{j},t_{k+1})}{\Delta x^{2}}]$

$\frac{\partial^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k+\frac{1}{2}})}\approx \frac{1}{2}[\frac{\partial^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k})}+\frac{\partial^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k+1})}]\approx\frac{1}{2}[\frac{u(x_{i},y_{j+1},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j-1},t_{k})}{\Delta y^{2}}+\frac{u(x_{i},y_{j+1},t_{k+1})-2u(x_{i},y_{j},t_{k+1})+u(x_{i},y_{j-1},t_{k+1})}{\Delta y^{2}}]$

将上面三式代入公式（4），用数值解代替精确解并忽略高阶项，可得数值格式为

$\frac{u^{k+1}_{i,j}-u^{k}_{i,j}}{\Delta t}-\frac{1}{2}(\frac{u^{k+1}_{i+1,j}-2u^{k+1}_{i,j}+u^{k+1}_{i-1,j}+u^{k}_{i+1,j}-2u^{k}_{i,j}+u^{k}_{i-1,j}}{\Delta x^{2}}+\frac{u^{k+1}_{i,j+1}-2u^{k+1}_{i,j}+u^{k+1}_{i,j-1}+u^{k}_{i,j+1}-2u^{k}_{i,j}+u^{k}_{i,j-1}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k+\frac{1}{2}}) \space\space(5)$

上式可整理为

$\left\{\begin{matrix} -\frac{r_{1}}{2}u^{k+1}_{i-1,j}+(1+r_{1}+r_{2})u^{k+1}_{i,j}-\frac{r_{1}}{2}u^{k+1}_{i+1,j}-\frac{r_{2}}{2}u^{k+1}_{i,j-1}-\frac{r_{2}}{2}u^{k+1}_{i,j+1}=\\ \frac{r_{1}}{2}u^{k}_{i-1,j}+(1-r_{1}-r_{2})u^{k}_{i,j}+\frac{r_{1}}{2}u^{k}_{i+1,j}+\frac{r_{2}}{2}u^{k}_{i,j-1}+\frac{r_{2}}{2}u^{k}_{i,j+1}+f(x_{i},y_{j},t_{k})\Delta t, \\1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,0\leqslant k\leqslant l-1, \\ u^{0}_{i,j}=\varphi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\\ u^{k}_{0,j}=g_{1}(y_{j},t_{k}),u^{k}_{m,j}=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u^{k}_{i,0}=g_{3}(x_{i},t_{k}),u^{k}_{i,n}=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.$

上式写成矩阵形式为

$\begin{pmatrix} -\frac{r_{2}}{2} & & & & \\ & -\frac{r_{2}}{2} & & 0 & \\ & & \vdots & & \\ & 0 & & -\frac{r_{2}}{2} & \\ & & & & -\frac{r_{2}}{2} \end{pmatrix}$ $\begin{pmatrix} u^{k+1}_{1,j-1}\\ u^{k+1}_{2,j-1}\\ \vdots\\ u^{k+1}_{m-2,j-1}\\ u^{k+1}_{m-1,j-1} \end{pmatrix}+$ $\begin{pmatrix} 1+r_{1}+r_{2} & -\frac{r_{1}}{2} & & & \\ -\frac{r_{1}}{2}& 1+r_{1}+r_{2} & -\frac{r_{1}}{2} & 0 & \\ & & \ddots & & \\ & 0 & -\frac{r_{1}}{2} & 1+r_{1}+r_{2} &-\frac{r_{1}}{2} \\ & & & -\frac{r_{1}}{2}& 1+r_{1}+r_{2} \end{pmatrix}$ $\begin{pmatrix} u^{k+1}_{1,j}\\ u^{k+1}_{2,j}\\ \vdots\\ u^{k+1}_{m-2,j}\\ u^{k+1}_{m-1,j} \end{pmatrix}+$ $\begin{pmatrix} -\frac{r_{2}}{2} & & & & \\ & -\frac{r_{2}}{2} & & 0 & \\ & & \vdots & & \\ & 0 & & -\frac{r_{2}}{2} & \\ & & & & -\frac{r_{2}}{2} \end{pmatrix}$ $\begin{pmatrix} u^{k+1}_{1,j+1}\\ u^{k+1}_{2,j+1}\\ \vdots\\ u^{k+1}_{m-2,j+1}\\ u^{k+1}_{m-1,j+1} \end{pmatrix}=$ $\begin{pmatrix} \frac{r_{1}}{2}u^{k}_{0,j}+(1-r_{1}-r_{2})u^{k}_{1,j}+\frac{r_{1}}{2}u^{k}_{2,j}+\frac{r_{2}}{2}u^{k}_{1,j-1}+\frac{r_{2}}{2}u^{k}_{1,j+1}+f(x_{1},y_{j},t_{k+\frac{1}{2}})\Delta t+\frac{r_{1}}{2}u^{k+1}_{0,j}\\ \frac{r_{1}}{2}u^{k}_{1,j}+(1-r_{1}-r_{2})u^{k}_{2,j}+\frac{r_{1}}{2}u^{k}_{3,j}+\frac{r_{2}}{2}u^{k}_{2,j-1}+\frac{r_{2}}{2}u^{k}_{2,j+1}+f(x_{2},y_{j},t_{k+\frac{1}{2}})\Delta t\\ \vdots \\ \frac{r_{1}}{2}u^{k}_{m-3,j}+(1-r_{1}-r_{2})u^{k}_{m-2,j}+\frac{r_{1}}{2}u^{k}_{m-1,j}+\frac{r_{2}}{2}u^{k}_{m-2,j-1}+\frac{r_{2}}{2}u^{k}_{m-2,j+1}+f(x_{m-2},y_{j},t_{k+\frac{1}{2}})\Delta t\\ \frac{r_{1}}{2}u^{k}_{m-2,j}+(1-r_{1}-r_{2})u^{k}_{m-1,j}+\frac{r_{1}}{2}u^{k}_{m,j}+\frac{r_{2}}{2}u^{k}_{m-1,j-1}+\frac{r_{2}}{2}u^{k}_{m-1,j+1}+f(x_{m-1},y_{j},t_{k+\frac{1}{2}})\Delta t+\frac{r_{1}}{2}u^{k+1}_{m,j} \end{pmatrix}$

上式为五对角矩阵的线性方程组，求解困难。

2.3 交替方向隐(ADI)格式

为了分别求解对角矩阵，我们需要用到一个二阶中心差分记号，即算子符合：

$\left\{\begin{matrix} \delta ^{2}_{x}u(x_{i},y_{j},t_{k})=u(x_{i-1},y_{j},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i+1},y_{j},t_{k})\\ \delta ^{2}_{y}u(x_{i},y_{j},t_{k})=u(x_{i},y_{j-1},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j+1},t_{k})\\ \delta ^{2}_{t}u(x_{i},y_{j},t_{k-1})=u(x_{i},y_{j},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j},t_{k+1}) \end{matrix}\right.$

及其离散形式：

$\left\{\begin{matrix} \delta^{2}_{x}u^{k}_{ij}=u^{k}_{i-1,j}-2u^{k}_{ij}+u^{k}_{i+1,j},\\ \delta^{2}_{y}u^{k}_{ij}=u^{k}_{i,j-1}-2u^{k}_{ij}+u^{k}_{i,j+1},\\ \delta^{2}_{t}u^{k}_{ij}=u^{k-1}_{ij}-2u^{k}_{ij}+u^{k+1}_{ij} \end{matrix}\right.$

利用这些记号，可以将公式（5）写成

$\frac{u^{k+1}_{i,j}-u^{k}_{i,j}}{\Delta t}-\frac{1}{2}(\frac{\delta^{2}_{x}u^{k+1}_{i,j}+\delta^{2}_{x}u^{k}_{i,j}}{\Delta x^{2}}+\frac{\delta^{2}_{y}u^{k+1}_{i,j}+\delta^{2}_{y}u^{k}_{i,j}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k+\frac{1}{2}})$

整理可得

$u^{k+1}_{i,j}-u^{k}_{i,j}-\frac{1}{2}(r_{1}(\delta^{2}_{x}u^{k+1}_{i,j}+\delta^{2}_{x}u^{k}_{i,j})+r_{2}(\delta^{2}_{y}u^{k+1}_{i,j}+\delta^{2}_{y}u^{k}_{i,j}))=f(x_{i},y_{j},t_{k+\frac{1}{2}})\Delta t$

再整理，可得

$(1-\frac{r_{1}}{2}\delta^{2}_{x}-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=(1+\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+f(x_{i},y_{j},t_{k+\frac{1}{2}})\Delta t \space\space(6)$

在公式（6）左右两端分别添加辅助项 $\frac{r_{1}r_{2}}{4}\delta^{2}_{x}\delta^{2}_{y}u^{k+1}_{i,j}$ 和 $\frac{r_{1}r_{2}}{4}\delta^{2}_{x}\delta^{2}_{y}u^{k}_{i,j}$ 。于是Crank-Nicolson格式可以修正为：

$(1-\frac{r_{1}}{2}\delta^{2}_{x}-\frac{r_{2}}{2}\delta^{2}_{y}+\frac{r_{1}r_{2}}{4}\delta^{2}_{x}\delta^{2}_{y})u^{k+1}_{i,j}=(1+\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}+\frac{r_{1}r_{2}}{4}\delta^{2}_{x}\delta^{2}_{y})u^{k}_{i,j}+f(x_{i},y_{j},t_{k+\frac{1}{2}})\Delta t$

即

$(1-\frac{r_{1}}{2}\delta^{2}_{x})(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=(1+\frac{r_{1}}{2}\delta^{2}_{x})(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+f^{k+\frac{1}{2}}_{ij}\Delta t$

其中， $f^{k+\frac{1}{2}}_{i,j}=f(x_{i},y_{j},t_{k+\frac{1}{2}})$ 。为使分解更彻底，上式可以写成

$(1-\frac{r_{1}}{2}\delta^{2}_{x})(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=(1+\frac{r_{1}}{2}\delta^{2}_{x})(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+\frac{1}{2}(1-\frac{r_{1}}{2}\delta^{2}_{x}+1+\frac{r_{1}}{2}\delta^{2}_{x})f^{k+\frac{1}{2}}_{i,j}\Delta t$

进一步分解得

$(1-\frac{r_{1}}{2}\delta^{2}_{x})[(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}-\frac{1}{2}f^{k+\frac{1}{2}}_{i,j}\Delta t]=(1+\frac{r_{1}}{2}\delta^{2}_{x})[(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+\frac{1}{2}f^{k+\frac{1}{2}}_{i,j}\Delta t] \space\space\space\space (7)$

引入中间变量 $V_{i,j}$ ，使公式（7）中最右端方括号内的元素满足：

$(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+\frac{1}{2}f^{k+\frac{1}{2}}_{i,j}\Delta t =(1-\frac{r_{1}}{2}\delta^{2}_{x})V_{i,j} \space\space(8)$

再由于差分算子具有可交换性，由公式（7）有：

$(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}-\frac{1}{2}f^{k+\frac{1}{2}}_{i,j}\Delta t =(1+\frac{r_{1}}{2}\delta^{2}_{x})V_{i,j} \space\space(9)$

2.3.1 Peaceman-Rachford格式

组合公式（8）和公式（9），可得Peaceman-Rachford格式：

$\left\{\begin{matrix} (1-\frac{r_{1}}{2}\delta^{2}_{x})V_{i,j}=(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+\frac{1}{2}f^{k+\frac{1}{2}}_{i,j}\Delta t \\ (1+\frac{r_{1}}{2}\delta^{2}_{x})V_{i,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}-\frac{1}{2}f^{k+\frac{1}{2}}_{i,j}\Delta t \end{matrix}\right.$ $(10)$

将公式（10）中两式相减，整理后可得：

$V_{i,j}=\frac{1}{2}(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+\frac{1}{2}(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}$

从而有：

$\left\{\begin{matrix} V_{0,j}=\frac{1}{2}(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{0,j}+\frac{1}{2}(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{0,j}\\ V_{m,j}=\frac{1}{2}(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{m,j}+\frac{1}{2}(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{m,j} \end{matrix}\right.$ $(11)$

公式（10）和（11）即为完整的Peaceman-Rachford格式。整个求解过程从k=0第0层开始逐层求解。这种方式是关于x方向和y方向交替求解的隐式方法，称为ADI方法。将其写成差分形式：

$\left\{\begin{matrix} -\frac{r_{1}}{2}V_{i-1,j}+(1+r_{1})V_{i,j}-\frac{r_{1}}{2}V_{i+1,j}=\frac{r_{2}}{2}u^{k}_{i,j-1}+(1-r_{2})u^{k}_{i,j}+\frac{r_{2}}{2}u^{k}_{i,j+1}+\frac{1}{2}f^{k+\frac{1}{2}}_{i,j}\Delta t\\ -\frac{r_{2}}{2}u^{k+1}_{i,j-1}+(1+r_{2})u^{k+1}_{i,j}-\frac{r_{2}}{2}u^{k+1}_{i,j+1}=\frac{r_{1}}{2}(V_{i-1,j}+V_{i+1,j})+(1-r_{1})V_{i,j}+\frac{1}{2}f^{k+\frac{1}{2}}_{i,j}\Delta t \end{matrix}\right.$

其中， $1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,0\leqslant k\leqslant l-1,$ 且

$\left\{\begin{matrix} V_{0,j}=\frac{1-r_{2}}{2}u^{k}_{0,j}+\frac{1+r_{2}}{2}u^{k+1}_{0,j}+\frac{r_{2}}{4}(u^{k}_{0,j-1}+u^{k}_{0,j+1}-u^{k+1}_{0,j-1}-u^{k+1}_{0,j+1})\\ V_{m,j}=\frac{1-r_{2}}{2}u^{k}_{m,j}+\frac{1+r_{2}}{2}u^{k+1}_{m,j}+\frac{r_{2}}{4}(u^{k}_{m,j-1}+u^{k}_{m,j+1}-u^{k+1}_{m,j-1}-u^{k+1}_{m,j+1}) \end{matrix}\right.$

2.3.2 D’Yakonov格式

公式 $(1-\frac{r_{1}}{2}\delta^{2}_{x})(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=(1+\frac{r_{1}}{2}\delta^{2}_{x})(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+f^{k+\frac{1}{2}}_{ij}\Delta t$ 还可以分解为：

$\left\{\begin{matrix} (1-\frac{r_{1}}{2}\delta^{2}_{x})V_{i,j}=(1+\frac{r_{1}}{2}\delta^{2}_{x})(1+\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j}+f^{k+\frac{1}{2}}_{i,j}\Delta t\\ (1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=V_{i,j} \end{matrix}\right.$

$V_{0,j},V_{m,j}$ 可以直接从上式的第2式中获得：

$\left\{\begin{matrix} V_{0,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{0,j} \\ V_{m,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{m,j} \end{matrix}\right.$

将其写成差分形式：

$\left\{\begin{matrix} -\frac{r_{1}}{2}V_{i-1,j}+(1+r_{1})V_{i,j}-\frac{r_{1}}{2}V_{i+1,j}=\frac{r_{1}(1-r_{2})}{2}(u^{k}_{i-1,j}+u^{k}_{i+1,j}) +\frac{r_{2}(1-r_{1})}{2}(u^{k}_{i,j-1}+u^{k}_{i,j+1}) +\\(1-r_{1})(1-r_{2})u^{k}_{i,j}+f^{k+\frac{1}{2}}_{i,j}\Delta t +\frac{r_{1}r_{2}}{4}(u^{k}_{i-1,j-1}+u^{k}_{i+1,j-1}+u^{k}_{i-1,j+1}+u^{k}_{i+1,j+1}) \\ -\frac{r_{2}}{2}u^{k+1}_{i,j-1}+(1+r_{2})u^{k+1}_{i,j}-\frac{r_{2}}{2}u^{k+1}_{i,j+1}=V_{i,j} \end{matrix}\right.$

其中， $1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,0\leqslant k\leqslant l-1,$ 且

$\left\{\begin{matrix} V_{0,j}=-\frac{r_{2}}{2}(u^{k+1}_{0,j-1}+u^{k+1}_{0,j+1})+(1+r_{2})u^{k+1}_{0,j} \\ V_{m,j}=-\frac{r_{2}}{2}(u^{k+1}_{m,j-1}+u^{k+1}_{m,j+1})+(1+r_{2})u^{k+1}_{m,j} \end{matrix}\right.$

2.3.3 Douglas格式

$(1-\frac{r_{1}}{2}\delta^{2}_{x})(1-\frac{r_{2}}{2}\delta^{2}_{y})(u^{k+1}_{i,j}-u^{k}_{i,j})=(r_{1}\delta^{2}_{x}+r_{2}\delta^{2}_{y})u^{k}_{i,j}+f^{k+\frac{1}{2}}_{i,j})\Delta t$

借助中间变量 $V_{i,j}$ 可以分解为下式：

$\left\{\begin{matrix} (1-\frac{r_{1}}{2}\delta^{2}_{x})V_{i,j}=(r_{1}\delta^{2}_{x}+r_{2}\delta^{2}_{y})u^{k}_{i,j}+f^{k+\frac{1}{2}}_{i,j}\Delta t \\ (1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=V_{i,j}+(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{i,j} \end{matrix}\right.$

$V_{0,j},V_{m,j}$ 可以直接从上式的第2式中获得：

$\left\{\begin{matrix} V_{0,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{0,j} -(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{0,j} \\ V_{m,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{m,j}-(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k}_{m,j} \end{matrix}\right.$

将其写成差分形式：

$\left\{\begin{matrix} -\frac{r_{1}}{2}V_{i-1,j}+(1+r_{1})V_{i,j}-\frac{r_{1}}{2}V_{i+1,j}=r_{1}(u^{k}_{i-1,j}-2u^{k}_{i,j}+u^{k}_{i+1,j}) +r_{2}(u^{k}_{i,j-1}-2u^{k}_{i,j}+u^{k}_{i,j+1}) +f^{k+\frac{1}{2}}_{i,j}\Delta t \\ -\frac{r_{2}}{2}u^{k+1}_{i,j-1}+(1+r_{2})u^{k+1}_{i,j}-\frac{r_{2}}{2}u^{k+1}_{i,j+1}=V_{i,j} -\frac{r_{2}}{2}(u^{k}_{i,j-1}+u^{k}_{i,j+1})+(1+r_{2})u^{k}_{i,j} \end{matrix}\right.$

其中， $1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,0\leqslant k\leqslant l-1,$ 且

$\left\{\begin{matrix} V_{0,j}= -\frac{r_{2}}{2}(u^{k+1}_{0,j-1}+u^{k+1}_{0,j+1}-u^{k}_{0,j-1}-u^{k}_{0,j+1})+(1+r_{2})(u^{k+1}_{0,j} - u^{k}_{0,j}) \\ V_{m,j}=-\frac{r_{2}}{2}(u^{k+1}_{m,j-1}+u^{k+1}_{m,j+1}-u^{k}_{m,j-1}-u^{k}_{m,j+1})+(1+r_{2})(u^{k+1}_{m,j} - u^{k}_{m,j}) \end{matrix}\right.$

三、算例实现

二维抛物型方程初边值问题：

$\left\{\begin{matrix} \frac{\partial u}{\partial t}-(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=-\frac{3}{2}e^{(x+y)/2-t},(x,y)\in [0,1]\times[0,1],0<t\leqslant 1,\\ u(x,y,0)=e^{(x+y)/2},(x,y)\in[0,1]\times[0,1],\\ u(0,y,t)=e^{y/2-t},u(1,y,t)=e^{(1+y)/2-t},0\leqslant y\leqslant 1,0<t\leqslant 1,\\ u(x,0,t)=e^{x/2-t},u(x,1,t)=e^{(1+x)/2-t},0\leqslant x\leqslant 1,0<t\leqslant 1 \end{matrix}\right.$

已知精确解为 $u(x,y,t)=e^{(x+y)/2-t}$ 。分别计算两种不同步长条件下的数值解，并输出18个节点处的数值解和误差。【18个点分别为：（0.25i，0.25j，0.50）及（0.25i，0.25j，1.00），其中i，j=1,2,3.】

条件1： $\Delta x =\frac{1}{60},\Delta y =\frac{1}{40},\Delta t =\frac{1}{20}$

条件2： $\Delta x =\frac{1}{120},\Delta y =\frac{1}{80},\Delta t =\frac{1}{40}$

3.1 Peaceman-Rachford格式

代码如下：

#include <cmath>
#include<stdlib.h>
#include<stdio.h>

int main(int argc, char* argv[])
{
        int i, j, k, m, n, L, gap_i, gap_j;
        double a, b, T, r1, r2, dx, dy, dt;
        double  *x, *y, *t, ***u, **v;
        double *a1, *b1, *c1, *d1, *a2, *b2, *c2, *d2, tmid, *ans, temp;
        double f(double x, double y, double t);
        double phi(double x, double y);
        double g1(double y, double t);
        double g2(double y, double t);
        double g3(double x, double t);
        double g4(double x, double t);
        double *chase_algorithm(double *a, double *b, double *c, double *d, int n);
        double exact(double x, double y, double t);

        a=1.0;
        b=1.0;
        T=1.0;
        m=60;
        n=40;
        L=20;
        dx=a/m;
        dy=b/n;
        dt=T/L;
        r1=dt/(dx*dx);
        r2=dt/(dy*dy);
        printf("m=%d, n=%d, L=%d\n", m, n, L);
        printf("r1=%.4f,   r2=%.4f\n", r1, r2);

        x=(double *)malloc(sizeof(double)*(m+1));
        for(i=0;i<=m;i++)
                x[i]=i*dx;

        y=(double *)malloc(sizeof(double)*(n+1));
        for(j=0;j<=n;j++)
                y[j]=j*dy;

        t=(double *)malloc(sizeof(double)*(L+1));
        for(k=0;k<=L;k++)
                t[k]=k*dt;

        u=(double ***)malloc(sizeof(double *)*((m+1)*(n+1)*(L+1)));
        for(i=0;i<=m;i++)
                u[i]=(double **)malloc(sizeof(double *)*((n+1)*(L+1)));
        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                        u[i][j]=(double *)malloc(sizeof(double)*(L+1));
        }
                                                         
        v=(double **)malloc(sizeof(double *)*(m+1));
        for(i=0;i<=m;i++)
                v[i]=(double *)malloc(sizeof(double)*(n+1));

        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                {
                        u[i][j][0]=phi(x[i], y[j]);
                }
        }

        for(k=1;k<=L;k++)
        {
                for(j=0;j<=n;j++)
                {
                         u[0][j][k]=g1(y[j], t[k]);
                         u[m][j][k]=g2(y[j], t[k]);
                }
                for(i=1;i<=m-1;i++)
                {
                         u[i][0][k]=g3(x[i], t[k]);
                         u[i][n][k]=g4(x[i], t[k]);
                }
        }

        a1=(double *)malloc(sizeof(double)*(m-1));
        b1=(double *)malloc(sizeof(double)*(m-1));
        c1=(double *)malloc(sizeof(double)*(m-1));
        d1=(double *)malloc(sizeof(double)*(m-1));
        for(i=0;i<m-1;i++)
        {
                a1[i]=-r1/2.0;
                b1[i]=1.0+r1;
                c1[i]=a1[i];
        }

        a2=(double *)malloc(sizeof(double)*(n-1));
        b2=(double *)malloc(sizeof(double)*(n-1));
        c2=(double *)malloc(sizeof(double)*(n-1));
        d2=(double *)malloc(sizeof(double)*(n-1));
        for(j=0;j<n-1;j++)
        {
                a2[j]=-r2/2.0;
                b2[j]=1.0+r2;
                c2[j]=a2[j];
        }

        for(k=0;k<L;k++)
        {
                tmid=(t[k]+t[k+1])/2.0;
                for(j=1;j<=n-1;j++)
                {
                        for(i=1;i<=m-1;i++)
                                  d1[i-1]=r2*(u[i][j-1][k]+u[i][j+1][k])/2.0+(1.0-r2)*u[i][j][k]+f(x[i], y[j], tmid)*dt/2.0;
                        v[0][j]=(1-r2)*u[0][j][k]/2.0+(1+r2)*u[0][j][k+1]/2.0+r2*(u[0][j-1][k]+u[0][j+1][k]-u[0][j-1][k+1]-u[0][j+1][k+1])/4.0;
                        v[m][j]=(1-r2)*u[m][j][k]/2.0+(1+r2)*u[m][j][k+1]/2.0+r2*(u[m][j-1][k]+u[m][j+1][k]-u[m][j-1][k+1]-u[m][j+1][k+1])/4.0;
                        d1[0]=d1[0]+r1*v[0][j]/2.0;
                        d1[m-2]=d1[m-2]+r1*v[m][j]/2.0;
                        ans=chase_algorithm(a1, b1, c1, d1, m-1);

                        for(i=1;i<=m-1;i++)
                                v[i][j]=ans[i-1];
                        free(ans);
                }

                for(i=1;i<=m-1;i++)
                {
                        for(j=1;j<=n-1;j++)
                                d2[j-1]=r1*(v[i-1][j]+v[i+1][j])/2.0+(1.0-r1)*v[i][j]+f(x[i], y[j], tmid)*dt/2.0;
                        d2[0]=d2[0]+r2*u[i][0][k+1]/2.0;
                        d2[n-2]=d2[n-2]+r2*u[i][n][k+1]/2.0;
                        ans=chase_algorithm(a2, b2, c2, d2, n-1);

                        for(j=1;j<=n-1;j++)
                                u[i][j][k+1]=ans[j-1];
                        free(ans);
                }
        }

        gap_i=m/4;
        gap_j=n/4;

        for(i=gap_i;i<=m-1;i=i+gap_i)
        {
                for(j=gap_j;j<=n-1;j=j+gap_j)
                {
                        temp=fabs(exact(x[i], y[j], T/2.0)-u[i][j][L/2]);
                        printf("(%.2f, %.2f, 0.50)     y=%f,    err=%.4e\n", x[i], y[j], u[i][j][L/2], temp);
                }
        }

        printf("\n");
        printf("\n");

        for(i=gap_i;i<=m-1;i=i+gap_i)
        {
                for(j=gap_j;j<=n-1;j=j+gap_j)
                {
                        temp=fabs(exact(x[i], y[j], T)-u[i][j][L]);
                        printf("(%.2f, %.2f, 1.00)     y=%f,    err=%.4e\n", x[i], y[j], u[i][j][L], temp);
                }
        }

        free(x); free(y); free(t);
        free(a1); free(b1); free(c1); free(d1);
        free(a2); free(b2); free(c2); free(d2);

        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                      free(u[i][j]);
        }
        for(i=0;i<=m;i++)
        {
                free(u[i]); free(v[i]);
        }

        free(u); free(v);

        return 0;
}



double f(double x, double y, double t)
{
        return -3.0*exp((x+y)/2.0-t)/2.0;
}
double phi(double x, double y)
{
        return exp((x+y)/2.0);
}
double g1(double y, double t)
{
        return exp(y/2.0-t);
}
double g2(double y, double t)
{
        return exp((1.0+y)/2.0-t);
}
double g3(double x, double t)
{
        return exp(x/2.0-t);
}
double g4(double x, double t)
{
        return exp((1.0+x)/2.0-t);
}
double * chase_algorithm(double *a, double *b, double *c, double *d, int n)
{
        int i;
        double *ans, *g, *w, p;

        ans=(double *)malloc(sizeof(double)*n);
        g=(double *)malloc(sizeof(double)*n);
        w=(double *)malloc(sizeof(double)*n);
        g[0]=d[0]/b[0];
        w[0]=c[0]/b[0];

        for(i=1;i<n;i++)
        {
                p=b[i]-a[i]*w[i-1];
                g[i]=(d[i]-a[i]*g[i-1])/p;
                w[i]=c[i]/p;
        }

        ans[n-1]=g[n-1];
        i=n-2;
        do
        {
                ans[i]=g[i]-w[i]*ans[i+1];
                i=i-1;
        }while(i>=0);

        free(g);free(w);
        return ans;
}
double exact(double x, double y, double t)
{
        return exp((x+y)/2.0-t);
}

当m=60，n=40，L=20时，结果如下：

m=60, n=40, L=20
r1=180.0000,   r2=80.0000
(0.25, 0.25, 0.50)     y=0.778814,    err=1.2876e-05
(0.25, 0.50, 0.50)     y=0.882514,    err=1.7519e-05
(0.25, 0.75, 0.50)     y=1.000015,    err=1.4611e-05
(0.50, 0.25, 0.50)     y=0.882514,    err=1.7519e-05
(0.50, 0.50, 0.50)     y=1.000024,    err=2.3999e-05
(0.50, 0.75, 0.50)     y=1.133168,    err=1.9757e-05
(0.75, 0.25, 0.50)     y=1.000015,    err=1.4611e-05
(0.75, 0.50, 0.50)     y=1.133168,    err=1.9757e-05
(0.75, 0.75, 0.50)     y=1.284042,    err=1.6707e-05


(0.25, 0.25, 1.00)     y=0.472374,    err=7.8105e-06
(0.25, 0.50, 1.00)     y=0.535272,    err=1.0627e-05
(0.25, 0.75, 1.00)     y=0.606540,    err=8.8626e-06
(0.50, 0.25, 1.00)     y=0.535272,    err=1.0627e-05
(0.50, 0.50, 1.00)     y=0.606545,    err=1.4557e-05
(0.50, 0.75, 1.00)     y=0.687301,    err=1.1984e-05
(0.75, 0.25, 1.00)     y=0.606540,    err=8.8626e-06
(0.75, 0.50, 1.00)     y=0.687301,    err=1.1984e-05
(0.75, 0.75, 1.00)     y=0.778811,    err=1.0134e-05

当m=120，n=80，L=40时，结果如下：

m=120, n=80, L=40
r1=360.0000,   r2=160.0000
(0.25, 0.25, 0.50)     y=0.778804,    err=3.2127e-06
(0.25, 0.50, 0.50)     y=0.882501,    err=4.3716e-06
(0.25, 0.75, 0.50)     y=1.000004,    err=3.6449e-06
(0.50, 0.25, 0.50)     y=0.882501,    err=4.3717e-06
(0.50, 0.50, 0.50)     y=1.000006,    err=5.9884e-06
(0.50, 0.75, 0.50)     y=1.133153,    err=4.9297e-06
(0.75, 0.25, 0.50)     y=1.000004,    err=3.6449e-06
(0.75, 0.50, 0.50)     y=1.133153,    err=4.9297e-06
(0.75, 0.75, 0.50)     y=1.284030,    err=4.1668e-06


(0.25, 0.25, 1.00)     y=0.472369,    err=1.9488e-06
(0.25, 0.50, 1.00)     y=0.535264,    err=2.6518e-06
(0.25, 0.75, 1.00)     y=0.606533,    err=2.2109e-06
(0.50, 0.25, 1.00)     y=0.535264,    err=2.6518e-06
(0.50, 0.50, 1.00)     y=0.606534,    err=3.6325e-06
(0.50, 0.75, 1.00)     y=0.687292,    err=2.9903e-06
(0.75, 0.25, 1.00)     y=0.606533,    err=2.2109e-06
(0.75, 0.50, 1.00)     y=0.687292,    err=2.9902e-06
(0.75, 0.75, 1.00)     y=0.778803,    err=2.5275e-06

3.2 D’Yakonov格式

代码如下：

#include <cmath>
#include <stdlib.h>
#include <stdio.h>


int main(int argc, char* argv[])
{
        int i, j, k, m, n, L, gap_i, gap_j;
        double a, b, T, r1, r2, dx, dy, dt;
        double *x, *y, *t, ***u, **v;
        double *a1, *b1, *c1, *d1, *a2, *b2, *c2, *d2, tmid, *ans, temp;
        double f(double x, double y, double t);
        double phi(double x, double y);
        double g1(double y, double t);
        double g2(double y, double t);
        double g3(double x, double t);
        double g4(double x, double t);
        double *chase_algorithm(double *a, double *b, double *c, double *d, int n);
        double exact(double x, double y, double t);

        a=1.0;
        b=1.0;
        T=1.0;
        m=60;
        n=40;
        L=20;
        dx=a/m;
        dy=b/n;
        dt=T/L;
        r1=dt/(dx*dx);
        r2=dt/(dy*dy);
        printf("m=%d, n=%d, L=%d\n", m, n, L);
        printf("r1=%.4f,   r2=%.4f\n", r1, r2);

        x=(double *)malloc(sizeof(double)*(m+1));
        for(i=0;i<=m;i++)
                x[i]=i*dx;

        y=(double *)malloc(sizeof(double)*(n+1));
        for(j=0;j<=n;j++)
                y[j]=j*dy;
        t=(double *)malloc(sizeof(double)*(L+1));
        for(k=0;k<=L;k++)
                t[k]=k*dt;

        u=(double ***)malloc(sizeof(double *)*((m+1)*(n+1)*(L+1)));
        for(i=0;i<=m;i++)
                u[i]=(double **)malloc(sizeof(double *)*((n+1)*(L+1)));
        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                u[i][j]=(double *)malloc(sizeof(double)*(L+1));
        }
        
        v=(double **)malloc(sizeof(double *)*(m+1));
        for(i=0;i<=m;i++)
                v[i]=(double *)malloc(sizeof(double)*(n+1));

        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                {
                        u[i][j][0]=phi(x[i], y[j]);
                }
        }

        for(k=1;k<=L;k++)
        {
                for(j=0;j<=n;j++)
                {
                        u[0][j][k]=g1(y[j], t[k]);
                        u[m][j][k]=g2(y[j], t[k]);
                }
                for(i=1;i<=m-1;i++)
                {
                        u[i][0][k]=g3(x[i], t[k]);
                        u[i][n][k]=g4(x[i], t[k]);
                }
        }

        a1=(double *)malloc(sizeof(double)*(m-1));
        b1=(double *)malloc(sizeof(double)*(m-1));
        c1=(double *)malloc(sizeof(double)*(m-1));
        d1=(double *)malloc(sizeof(double)*(m-1));
        for(i=0;i<m-1;i++)
        {
                a1[i]=-r1/2.0;
                b1[i]=1.0+r1;
                c1[i]=a1[i];
        }

        a2=(double *)malloc(sizeof(double)*(n-1));
        b2=(double *)malloc(sizeof(double)*(n-1));
        c2=(double *)malloc(sizeof(double)*(n-1));
        d2=(double *)malloc(sizeof(double)*(n-1));
        for(j=0;j<n-1;j++)
        {
                a2[j]=-r2/2.0;
                b2[j]=1.0+r2;
                c2[j]=a2[j];
        }

        for(k=0;k<L;k++)
        {
                tmid=(t[k]+t[k+1])/2.0;
                for(j=1;j<=n-1;j++)
                {
                        for(i=1;i<=m-1;i++)
                        {
                                temp=r2*(1-r1)*(u[i][j-1][k]+u[i][j+1][k])/2.0+r1*(1-r2)*(u[i-1][j][k]+u[i+1][j][k])/2.0+(1-r1)*(1-r2)*u[i][j][k];
                                d1[i-1]=temp+f(x[i], y[j], tmid)*dt+r1*r2*(u[i-1][j-1][k]+u[i+1][j-1][k]+u[i-1][j+1][k]+u[i+1][j+1][k])/4.0;
                        }
                        v[0][j]=(1+r2)*u[0][j][k+1]-r2*(u[0][j-1][k+1]+u[0][j+1][k+1])/2.0;
                        v[m][j]=(1+r2)*u[m][j][k+1]-r2*(u[m][j-1][k+1]+u[m][j+1][k+1])/2.0;
                        d1[0]=d1[0]+r1*v[0][j]/2.0;
                        d1[m-2]=d1[m-2]+r1*v[m][j]/2.0;
                        ans=chase_algorithm(a1, b1, c1, d1, m-1);

                        for(i=1;i<=m-1;i++)
                                v[i][j]=ans[i-1];
                        free(ans);
                }

                for(i=1;i<=m-1;i++)
                {
                        for(j=1;j<=n-1;j++)
                                d2[j-1]=v[i][j];
                        d2[0]=d2[0]+r2*u[i][0][k+1]/2.0;
                        d2[n-2]=d2[n-2]+r2*u[i][n][k+1]/2.0;
                        ans=chase_algorithm(a2, b2, c2, d2, n-1);

                        for(j=1;j<=n-1;j++)
                                u[i][j][k+1]=ans[j-1];

                        free(ans);
                }
        }

        gap_i=m/4;
        gap_j=n/4;

        for(i=gap_i;i<=m-1;i=i+gap_i)
        {
                for(j=gap_j;j<=n-1;j=j+gap_j)
                {
                        temp=fabs(exact(x[i], y[j], T/2.0)-u[i][j][L/2]);
                        printf("(%.2f, %.2f, 0.50)     y=%f,    err=%.4e\n", x[i], y[j], u[i][j][L/2], temp);
                }
        }

        printf("\n");
        printf("\n");

        for(i=gap_i;i<=m-1;i=i+gap_i)
        {
                for(j=gap_j;j<=n-1;j=j+gap_j)
                {
                        temp=fabs(exact(x[i], y[j], T)-u[i][j][L]);
                        printf("(%.2f, %.2f, 1.00)     y=%f,    err=%.4e\n", x[i], y[j], u[i][j][L], temp);
                }
        }

        free(x); free(y); free(t);
        free(a1); free(b1); free(c1); free(d1);
        free(a2); free(b2); free(c2); free(d2);
        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                        free(u[i][j]);
        }

        for(i=0;i<=m;i++)
        {
                free(u[i]);free(v[i]);
        }

        free(u); free(v);

        return 0;
}


double f(double x, double y, double t)
{
        return -3.0*exp((x+y)/2.0-t)/2.0;
}
double phi(double x, double y)
{
        return exp((x+y)/2.0);
}
double g1(double y, double t)
{
        return exp(y/2.0-t);
}
double g2(double y, double t)
{
        return exp((1.0+y)/2.0-t);
}
double g3(double x, double t)
{
        return exp(x/2.0-t);
}
double g4(double x, double t)
{
        return exp((1.0+x)/2.0-t);
}
double *chase_algorithm(double *a, double *b, double *c, double *d, int n)
{
        int i;
        double *ans, *g, *w, p;

        ans=(double *)malloc(sizeof(double)*n);
        g=(double *)malloc(sizeof(double)*n);
        w=(double *)malloc(sizeof(double)*n);
        g[0]=d[0]/b[0];
        w[0]=c[0]/b[0];

        for(i=1;i<n;i++)
        {
                p=b[i]-a[i]*w[i-1];
                g[i]=(d[i]-a[i]*g[i-1])/p;
                w[i]=c[i]/p;
        }

        ans[n-1]=g[n-1];
        i=n-2;
        do
        {
                ans[i]=g[i]-w[i]*ans[i+1];
                i=i-1;
        }while(i>=0);

        free(g); free(w);

        return ans;
}
double exact(double x, double y, double t)
{
        return exp((x+y)/2.0-t);
}

当m=60，n=40，L=20时，结果如下：

m=60, n=40, L=20
r1=180.0000,   r2=80.0000
(0.25, 0.25, 0.50)     y=0.778814,    err=1.2876e-05
(0.25, 0.50, 0.50)     y=0.882514,    err=1.7519e-05
(0.25, 0.75, 0.50)     y=1.000015,    err=1.4611e-05
(0.50, 0.25, 0.50)     y=0.882514,    err=1.7519e-05
(0.50, 0.50, 0.50)     y=1.000024,    err=2.3999e-05
(0.50, 0.75, 0.50)     y=1.133168,    err=1.9757e-05
(0.75, 0.25, 0.50)     y=1.000015,    err=1.4611e-05
(0.75, 0.50, 0.50)     y=1.133168,    err=1.9757e-05
(0.75, 0.75, 0.50)     y=1.284042,    err=1.6707e-05


(0.25, 0.25, 1.00)     y=0.472374,    err=7.8105e-06
(0.25, 0.50, 1.00)     y=0.535272,    err=1.0627e-05
(0.25, 0.75, 1.00)     y=0.606540,    err=8.8626e-06
(0.50, 0.25, 1.00)     y=0.535272,    err=1.0627e-05
(0.50, 0.50, 1.00)     y=0.606545,    err=1.4557e-05
(0.50, 0.75, 1.00)     y=0.687301,    err=1.1984e-05
(0.75, 0.25, 1.00)     y=0.606540,    err=8.8626e-06
(0.75, 0.50, 1.00)     y=0.687301,    err=1.1984e-05
(0.75, 0.75, 1.00)     y=0.778811,    err=1.0134e-05

当m=120，n=80，L=40时，结果如下：

m=120, n=80, L=40
r1=360.0000,   r2=160.0000
(0.25, 0.25, 0.50)     y=0.778804,    err=3.2127e-06
(0.25, 0.50, 0.50)     y=0.882501,    err=4.3716e-06
(0.25, 0.75, 0.50)     y=1.000004,    err=3.6449e-06
(0.50, 0.25, 0.50)     y=0.882501,    err=4.3717e-06
(0.50, 0.50, 0.50)     y=1.000006,    err=5.9884e-06
(0.50, 0.75, 0.50)     y=1.133153,    err=4.9297e-06
(0.75, 0.25, 0.50)     y=1.000004,    err=3.6449e-06
(0.75, 0.50, 0.50)     y=1.133153,    err=4.9297e-06
(0.75, 0.75, 0.50)     y=1.284030,    err=4.1668e-06


(0.25, 0.25, 1.00)     y=0.472369,    err=1.9488e-06
(0.25, 0.50, 1.00)     y=0.535264,    err=2.6518e-06
(0.25, 0.75, 1.00)     y=0.606533,    err=2.2109e-06
(0.50, 0.25, 1.00)     y=0.535264,    err=2.6518e-06
(0.50, 0.50, 1.00)     y=0.606534,    err=3.6325e-06
(0.50, 0.75, 1.00)     y=0.687292,    err=2.9903e-06
(0.75, 0.25, 1.00)     y=0.606533,    err=2.2109e-06
(0.75, 0.50, 1.00)     y=0.687292,    err=2.9902e-06
(0.75, 0.75, 1.00)     y=0.778803,    err=2.5275e-06

3.3 Douglas格式

代码如下：

#include <stdlib.h>
#include <stdio.h>

int main(int argc, char* argv[])
{
        int i, j, k, m, n, L, gap_i, gap_j;
        double a, b, T, r1, r2, dx, dy, dt;
        double *x, *y, *t, ***u, **v;
        double *a1, *b1, *c1, *d1, *a2, *b2, *c2, *d2, tmid, *ans, temp;
        double f(double x, double y, double t);
        double phi(double x, double y);
        double g1(double y, double t);
        double g2(double y, double t);
        double g3(double x, double t);
        double g4(double x, double t);
        double *chase_algorithm(double *a, double *b, double *c, double *d, int n);
        double exact(double x, double y, double t);

        a=1.0;
        b=1.0;
        T=1.0;
        m=60;
        n=40;
        L=20;
        dx=a/m;
        dy=b/n;
        dt=T/L;
        r1=dt/(dx*dx);
        r2=dt/(dy*dy);
        printf("m=%d, n=%d, L=%d\n", m, n, L);
        printf("r1=%.4f,   r2=%.4f\n", r1, r2);

        x=(double *)malloc(sizeof(double)*(m+1));
        for(i=0;i<=m;i++)
                x[i]=i*dx;

        y=(double *)malloc(sizeof(double)*(n+1));
        for(j=0;j<=n;j++)
                y[j]=j*dy;

        t=(double *)malloc(sizeof(double)*(L+1));
        for(k=0;k<=L;k++)
                t[k]=k*dt;

        u=(double ***)malloc(sizeof(double *)*((m+1)*(n+1)*(L+1)));
        for(i=0;i<=m;i++)
                u[i]=(double **)malloc(sizeof(double *)*((n+1)*(L+1)));
        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                        u[i][j]=(double *)malloc(sizeof(double)*(L+1));
        }

        v=(double **)malloc(sizeof(double *)*(m+1));
        for(i=0;i<=m;i++)
                v[i]=(double *)malloc(sizeof(double)*(n+1));

        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                {
                         u[i][j][0]=phi(x[i], y[j]);
                }
        }

        for(k=1;k<=L;k++)
        {
                for(j=0;j<=n;j++)
                {
                         u[0][j][k]=g1(y[j], t[k]);
                         u[m][j][k]=g2(y[j], t[k]);
                }
                for(i=1;i<=m-1;i++)
                {
                         u[i][0][k]=g3(x[i], t[k]);
                         u[i][n][k]=g4(x[i], t[k]);
                }
        }

        a1=(double *)malloc(sizeof(double)*(m-1));
        b1=(double *)malloc(sizeof(double)*(m-1));
        c1=(double *)malloc(sizeof(double)*(m-1));
        d1=(double *)malloc(sizeof(double)*(m-1));
        for(i=0;i<m-1;i++)
        {
                a1[i]=-r1/2.0;
                b1[i]=1.0+r1;
                c1[i]=a1[i];
        }

        a2=(double *)malloc(sizeof(double)*(n-1));
        b2=(double *)malloc(sizeof(double)*(n-1));
        c2=(double *)malloc(sizeof(double)*(n-1));
        d2=(double *)malloc(sizeof(double)*(n-1));
        for(j=0;j<n-1;j++)
        {
                a2[j]=-r2/2.0;
                b2[j]=1.0+r2;
                c2[j]=a2[j];
        }

        for(k=0;k<L;k++)
        {
                tmid=(t[k]+t[k+1])/2.0;
                for(j=1;j<=n-1;j++)
                {
                        for(i=1;i<=m-1;i++)
                        {
                                d1[i-1]=r1*(u[i-1][j][k]-2*u[i][j][k]+u[i+1][j][k])+r2*(u[i][j-1][k]-2*u[i][j][k]+u[i][j+1][k])+f(x[i], y[j], tmid)*dt;
                        }
                        v[0][j]=(1+r2)*(u[0][j][k+1]-u[0][j][k])-r2*(u[0][j+1][k+1]-u[0][j+1][k]+u[0][j-1][k+1]-u[0][j-1][k])/2.0;
                        v[m][j]=(1+r2)*(u[m][j][k+1]-u[m][j][k])-r2*(u[m][j+1][k+1]-u[m][j+1][k]+u[m][j-1][k+1]-u[m][j-1][k])/2.0;
                        d1[0]=d1[0]+r1*v[0][j]/2.0;
                        d1[m-2]=d1[m-2]+r1*v[m][j]/2.0;
                        ans=chase_algorithm(a1, b1, c1, d1, m-1);

                        for(i=1;i<=m-1;i++)
                                v[i][j]=ans[i-1];
                        free(ans);
                }

                for(i=1;i<=m-1;i++)
                {
                        for(j=1;j<=n-1;j++)
                                d2[j-1]=(1+r2)*u[i][j][k]-r2*(u[i][j-1][k]+u[i][j+1][k])/2.0+v[i][j];
                        d2[0]=d2[0]+r2*u[i][0][k+1]/2.0;
                        d2[n-2]=d2[n-2]+r2*u[i][n][k+1]/2.0;
                        ans=chase_algorithm(a2, b2, c2, d2, n-1);

                        for(j=1;j<=n-1;j++)
                                u[i][j][k+1]=ans[j-1];
                        free(ans);
                }
        }

        gap_i=m/4;
        gap_j=n/4;

        for(i=gap_i;i<=m-1;i=i+gap_i)
        {
                for(j=gap_j;j<=n-1;j=j+gap_j)
                {
                        temp=fabs(exact(x[i], y[j], T/2.0)-u[i][j][L/2]);
                        printf("(%.2f, %.2f, 0.50)     y=%f,    err=%.4e\n", x[i], y[j], u[i][j][L/2], temp);
                }
        }

        printf("\n");
        printf("\n");

        for(i=gap_i;i<=m-1;i=i+gap_i)
        {
               for(j=gap_j;j<=n-1;j=j+gap_j)
               {
                        temp=fabs(exact(x[i], y[j], T)-u[i][j][L]);
                        printf("(%.2f, %.2f, 1.00)     y=%f,    err=%.4e\n", x[i], y[j], u[i][j][L], temp);
               }
        }

        free(x); free(y); free(t);
        free(a1); free(b1); free(c1); free(d1);
        free(a2); free(b2); free(c2); free(d2);

        for(i=0;i<=m;i++)
        {
               for(j=0;j<=n;j++)
                       free(u[i][j]);
        }

        for(i=0;i<=m;i++)
        {
               free(u[i]); free(v[i]);
        }

        free(u); free(v);

        return 0;
}


double f(double x, double y, double t)
{
        return -3.0*exp((x+y)/2.0-t)/2.0;
}
double phi(double x, double y)
{
        return exp((x+y)/2.0);
}
double g1(double y, double t)
{
        return exp(y/2.0-t);
}
double g2(double y, double t)
{
        return exp((1.0+y)/2.0-t);
}
double g3(double x, double t)
{
        return exp(x/2.0-t);
}
double g4(double x, double t)
{
        return exp((1.0+x)/2.0-t);
}
double *chase_algorithm(double *a, double *b, double *c, double *d, int n)
{
        int i;
        double * ans, *g, *w, p;

        ans=(double *)malloc(sizeof(double)*n);
        g=(double *)malloc(sizeof(double)*n);
        w=(double *)malloc(sizeof(double)*n);

        g[0]=d[0]/b[0];
        w[0]=c[0]/b[0];

        for(i=1;i<n;i++)
        {
                p=b[i]-a[i]*w[i-1];
                g[i]=(d[i]-a[i]*g[i-1])/p;
                w[i]=c[i]/p;
        }

        ans[n-1]=g[n-1];
        i=n-2;
        do
        {
                ans[i]=g[i]-w[i]*ans[i+1];
                i=i-1;
        }while(i>=0);

        free(g); free(w);

        return ans;
}
double exact(double x, double y, double t)
{
        return exp((x+y)/2.0-t);
}

当m=60，n=40，L=20时，结果如下：

m=60, n=40, L=20
r1=180.0000,   r2=80.0000
(0.25, 0.25, 0.50)     y=0.778814,    err=1.2876e-05
(0.25, 0.50, 0.50)     y=0.882514,    err=1.7519e-05
(0.25, 0.75, 0.50)     y=1.000015,    err=1.4611e-05
(0.50, 0.25, 0.50)     y=0.882514,    err=1.7519e-05
(0.50, 0.50, 0.50)     y=1.000024,    err=2.3999e-05
(0.50, 0.75, 0.50)     y=1.133168,    err=1.9757e-05
(0.75, 0.25, 0.50)     y=1.000015,    err=1.4611e-05
(0.75, 0.50, 0.50)     y=1.133168,    err=1.9757e-05
(0.75, 0.75, 0.50)     y=1.284042,    err=1.6707e-05


(0.25, 0.25, 1.00)     y=0.472374,    err=7.8105e-06
(0.25, 0.50, 1.00)     y=0.535272,    err=1.0627e-05
(0.25, 0.75, 1.00)     y=0.606540,    err=8.8626e-06
(0.50, 0.25, 1.00)     y=0.535272,    err=1.0627e-05
(0.50, 0.50, 1.00)     y=0.606545,    err=1.4557e-05
(0.50, 0.75, 1.00)     y=0.687301,    err=1.1984e-05
(0.75, 0.25, 1.00)     y=0.606540,    err=8.8626e-06
(0.75, 0.50, 1.00)     y=0.687301,    err=1.1984e-05
(0.75, 0.75, 1.00)     y=0.778811,    err=1.0134e-05

当m=120，n=80，L=40时，结果如下：

m=120, n=80, L=40
r1=360.0000,   r2=160.0000
(0.25, 0.25, 0.50)     y=0.778804,    err=3.2127e-06
(0.25, 0.50, 0.50)     y=0.882501,    err=4.3716e-06
(0.25, 0.75, 0.50)     y=1.000004,    err=3.6449e-06
(0.50, 0.25, 0.50)     y=0.882501,    err=4.3717e-06
(0.50, 0.50, 0.50)     y=1.000006,    err=5.9884e-06
(0.50, 0.75, 0.50)     y=1.133153,    err=4.9297e-06
(0.75, 0.25, 0.50)     y=1.000004,    err=3.6449e-06
(0.75, 0.50, 0.50)     y=1.133153,    err=4.9297e-06
(0.75, 0.75, 0.50)     y=1.284030,    err=4.1668e-06


(0.25, 0.25, 1.00)     y=0.472369,    err=1.9488e-06
(0.25, 0.50, 1.00)     y=0.535264,    err=2.6518e-06
(0.25, 0.75, 1.00)     y=0.606533,    err=2.2109e-06
(0.50, 0.25, 1.00)     y=0.535264,    err=2.6518e-06
(0.50, 0.50, 1.00)     y=0.606534,    err=3.6325e-06
(0.50, 0.75, 1.00)     y=0.687292,    err=2.9903e-06
(0.75, 0.25, 1.00)     y=0.606533,    err=2.2109e-06
(0.75, 0.50, 1.00)     y=0.687292,    err=2.9902e-06
(0.75, 0.75, 1.00)     y=0.778803,    err=2.5275e-06