### 7.1 实验题目
用ADI法求解下列二维抛物方程初边值问题的解:
\[
\begin{cases}
\frac{\partial u}{\partial t}=\frac{1}{4^2}(u_{xx}+u_{yy}), (x,y)\in G=(0,1)\times(0,1),t > 0,\\
u(0,y,t)=u(1,y,t)=0, 0 < y < 1,t > 0,\\
u_y(x,0,t)=u_y(x,1,t)=0, 0 < x < 1,t > 0,\\
u(x,y,0)=\sin\pi x\cos\pi y.
\end{cases}
\]
精确解为:$u(x,y,t)=\sin\pi x\cos\pi y\exp(-\frac{\pi^2}{8}t)$。取空间步长$h = h_1 = h_2=\frac{1}{40}$,时间步长$\tau=\frac{1}{1600}$,网比$r = \frac{\tau}{h^2}=1$。计算到时间层$t = 1$,列出在节点$(x_j,y_k)$,$j,k = 1,2,3$的计算结果。并画出精确解曲面和数值解曲面。
### 7.2 实验原理
#### 第一步:网格剖分
设$x_j = jh(j = 0,1,\cdots,J)$,$y_k = kh(k = 0,1,\cdots,K)$,$t_n = n\tau(n = 0,1,\cdots,N)$差分解为$u^n_{j,k}$,则边值条件为:
\[
\begin{cases}
u^n_{0,k}=u^n_{J,k}=0, k = 0,\cdots,K\\
u^n_{j,0}=u^n_{j,1},u^n_{j,k - 1}=u^n_{j,k}, j = 0,1,\cdots,J.
\end{cases}
\]
初值条件为:$u^0_{j,k}=\sin\pi x_j\cos\pi y_k$。
#### 第二步:计算过程
从$n$到$n + 1$时,根据边值条件$u_{0,k}=u_{j,k}=0,k = 0,1,\cdots,K$,已经知道第0列和第$K$列数值全为0。
(i) 从$n\rightarrow n+\frac{1}{2}$,求$u^{n+\frac{1}{2}}_{j,k}$对$u_{xx}$向后差分,$u_{xx}$向前差分,构造出差分格式为:
\[
\begin{align*}
\frac{u^{n+\frac{1}{2}}_{j,k}-u^n_{j,k}}{\frac{\tau}{2}}&=\frac{1}{16}\left(\frac{u^{n+\frac{1}{2}}_{j + 1,k}-2u^{n+\frac{1}{2}}_{j,k}+u^{n+\frac{1}{2}}_{j - 1,k}}{h^2}+\frac{u^n_{j,k + 1}-2u^n_{j,k}+u^n_{j,k - 1}}{h^2}\right)\\
&=\frac{1}{16h^2}(\delta^2_xu^{n+\frac{1}{2}}_{j,k}+\delta^2_yu^n_{j,k})
\end{align*}
\]
从而得到:
\[
-\frac{1}{32}ru^{n+\frac{1}{2}}_{j + 1,k}+(1+\frac{1}{16}r)u^{n+\frac{1}{2}}_{j,k}-\frac{1}{32}ru^{n+\frac{1}{2}}_{j - 1,k}=\frac{1}{32}ru^n_{j,k + 1}+(1-\frac{1}{16}r)u^n_{j,k}+\frac{1}{32}ru^n_{j,k - 1}
\]
$j = 1,2,\cdots,J - 1$,$k = 1,2\cdots,K - 1$。
下面讨论$k = 1$时的情:
当$j = 1$,$-\frac{1}{32}ru^{n+\frac{1}{2}}_{2,1}+(1+\frac{1}{16}r)u^{n+\frac{1}{2}}_{1,1}-\frac{1}{32}ru^{n+\frac{1}{2}}_{0,1}=\frac{1}{32}ru^n_{1,2}+(1-\frac{1}{16}r)u^n_{1,1}+\frac{1}{32}ru^n_{1,0}$,
当$j = 2$,$-\frac{1}{32}ru^{n+\frac{1}{2}}_{3,1}+(1+\frac{1}{16}r)u^{n+\frac{1}{2}}_{2,1}-\frac{1}{32}ru^{n+\frac{1}{2}}_{1,1}=\frac{1}{32}ru^n_{2,2}+(1-\frac{1}{16}r)u^n_{2,1}+\frac{1}{32}ru^n_{2,0}$,
……
当$j = J - 2$,$-\frac{1}{32}ru^{n+\frac{1}{2}}_{J - 1,1}+(1+\frac{1}{16}r)u^{n+\frac{1}{2}}_{J - 2,1}-\frac{1}{32}ru^{n+\frac{1}{2}}_{J - 3,1}=\frac{1}{32}ru^n_{J - 2,2}+(1-\frac{1}{16}r)u^n_{J - 2,1}+\frac{1}{32}ru^n_{J - 2,0}$,
当$j = J - 1$,$-\frac{1}{32}ru^{n+\frac{1}{2}}_{J,1}+(1+\frac{1}{16}r)u^{n+\frac{1}{2}}_{J - 1,1}-\frac{1}{32}ru^{n+\frac{1}{2}}_{J - 2,1}=\frac{1}{32}ru^n_{J - 1,2}+(1-\frac{1}{16}r)u^n_{J - 1,1}+\frac{1}{32}ru^n_{J - 1,0}$。
再由边界条件得到$u^{n+\frac{1}{2}}_{0,1}=u^{n+\frac{1}{2}}_{J,1}=0$,则$k = 1$时,可以得到如下线性方程组:
\[
\begin{bmatrix}
1+\frac{1}{16}r&-\frac{1}{32}r&&&&\\
-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r&&&\\
&-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r&&\\
&&\ddots&\ddots&\ddots&\\
&&&-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r\\
&&&&-\frac{1}{32}r&1+\frac{1}{16}r
\end{bmatrix}
\begin{bmatrix}
u^{n+\frac{1}{2}}_{1,1}\\
u^{n+\frac{1}{2}}_{2,1}\\
u^{n+\frac{1}{2}}_{3,1}\\
\vdots\\
u^{n+\frac{1}{2}}_{J - 3,1}\\
u^{n+\frac{1}{2}}_{J - 2,1}\\
u^{n+\frac{1}{2}}_{J - 1,1}
\end{bmatrix}
=
\begin{bmatrix}
g_{1,1}\\
g_{2,1}\\
g_{3,1}\\
\vdots\\
g_{J - 3,1}\\
g_{J - 2,1}\\
g_{J - 1,1}
\end{bmatrix}
\]
则$k = 2$时,可以得到如下线性方程组:
\[
\begin{bmatrix}
1+\frac{1}{16}r&-\frac{1}{32}r&&&&\\
-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r&&&\\
&-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r&&\\
&&\ddots&\ddots&\ddots&\\
&&&-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r\\
&&&&-\frac{1}{32}r&1+\frac{1}{16}r
\end{bmatrix}
\begin{bmatrix}
u^{n+\frac{1}{2}}_{1,2}\\
u^{n+\frac{1}{2}}_{2,2}\\
u^{n+\frac{1}{2}}_{3,2}\\
\vdots\\
u^{n+\frac{1}{2}}_{J - 3,2}\\
u^{n+\frac{1}{2}}_{J - 2,2}\\
u^{n+\frac{1}{2}}_{J - 1,2}
\end{bmatrix}
=
\begin{bmatrix}
g_{1,2}\\
g_{2,2}\\
g_{3,2}\\
\vdots\\
g_{J - 3,2}\\
g_{J - 2,2}\\
g_{J - 1,2}
\end{bmatrix}
\]
依此类推:则$k = K - 1$时,可以得到如下线性方程组:
\[
\begin{bmatrix}
1+\frac{1}{16}r&-\frac{1}{32}r&&&&\\
-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r&&&\\
&-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r&&\\
&&\ddots&\ddots&\ddots&\\
&&&-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r\\
&&&&-\frac{1}{32}r&1+\frac{1}{16}r
\end{bmatrix}
\begin{bmatrix}
u^{n+\frac{1}{2}}_{1,K - 1}\\
u^{n+\frac{1}{2}}_{2,K - 1}\\
u^{n+\frac{1}{2}}_{3,K - 1}\\
\vdots\\
u^{n+\frac{1}{2}}_{J - 3,K - 1}\\
u^{n+\frac{1}{2}}_{J - 2,K - 1}\\
u^{n+\frac{1}{2}}_{J - 1,K - 1}
\end{bmatrix}
=
\begin{bmatrix}
g_{1,K - 1}\\
g_{2,K - 1}\\
g_{3,K - 1}\\
\vdots\\
g_{J - 3,K - 1}\\
g_{J - 2,K - 1}\\
g_{J - 1,K - 1}
\end{bmatrix}
\]
将上述方程组写成矩阵形式$AU^{n+\frac{1}{2}} = G$,其中
\[
A=\begin{bmatrix}
1+\frac{1}{16}r&-\frac{1}{32}r&&&&\\
-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r&&&\\
&-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r&&\\
&&\ddots&\ddots&\ddots&\\
&&&-\frac{1}{32}r&1+\frac{1}{16}r&-\frac{1}{32}r\\
&&&&-\frac{1}{32}r&1+\frac{1}{16}r
\end{bmatrix}
\]
\[
U^{n+\frac{1}{2}}=\begin{bmatrix}
u^{(n + 1/2)}_{1,1}&u^{(n + 1/2)}_{1,2}&\cdots&u^{(n + 1/2)}_{1,J - 1}\\
u^{(n + 1/2)}_{2,1}&u^{(n + 1/2)}_{2,2}&\cdots&u^{(n + 1/2)}_{2,J - 1}\\
u^{(n + 1/2)}_{3,1}&u^{(n + 1/2)}_{3,2}&\cdots&u^{(n + 1/2)}_{3,J - 1}\\
\vdots&\vdots&\ddots&\vdots\\
u^{(n + 1/2)}_{J - 3,1}&u^{(n + 1/2)}_{J - 3,2}&\cdots&u^{(n + 1/2)}_{J - 3,J - 1}\\
u^{(n + 1/2)}_{J - 2,1}&u^{(n + 1/2)}_{J - 2,2}&\cdots&u^{(n + 1/2)}_{J - 2,J - 1}\\
u^{(n + 1/2)}_{J - 1,1}&u^{(n + 1/2)}_{J - 1,2}&\cdots&u^{(n + 1/2)}_{J - 1,J - 1}
\end{bmatrix}
\]
\[
G=\begin{bmatrix}
g_{11}&g_{12}&\cdots&g_{1,J - 1}\\
g_{21}&g_{22}&\cdots&g_{2,J - 1}\\
g_{31}&g_{32}&\cdots&g_{3,J - 1}\\
\vdots&\vdots&\ddots&\vdots\\
g_{J - 3,1}&g_{J - 3,2}&\cdots&g_{J - 3,J - 1}\\
g_{J - 2,1}&g_{J - 2,2}&\cdots&g_{J - 2,J - 1}\\
g_{J - 1,1}&g_{J - 1,2}&\cdots&g_{J - 1,J - 1}
\end{bmatrix}
\]
又根据边值条件得:$u^{n+\frac{1}{2}}_{j,0}=u^{n+\frac{1}{2}}_{j,1},u^{n+\frac{1}{2}}_{j,K - 1}=u^{n+\frac{1}{2}}_{j,K},j = 0,1,\cdots,J$,解出第0行$u^{n+\frac{1}{2}}_{j,0}$和第$K$行$u^{n+\frac{1}{2}}_{j,K}(j = 1,2,\cdots,J - 1)$。而$j = 0$与$j = J$处的值由边值条件$u^{n+\frac{1}{2}}_{0,k}=u^{n+\frac{1}{2}}_{J,k}=0$得到。
(i) 第二步:从$n+\frac{1}{2}\rightarrow n + 1$,求$u^{n+\frac{1}{2}}_{j,k}$对$u_{xx}$向前差分,$u_{yy}$向后差分,构造出差分格式为:
\[
\begin{align*}
\frac{u^{n + 1}_{j,k}-u^n_{j,k}}{\frac{\tau}{2}}&=\frac{1}{16}\left(\frac{u^{n+\frac{1}{2}}_{j + 1,k}-2u^{n+\frac{1}{2}}_{j,k}+u^{n+\frac{1}{2}}_{j - 1,k}}{h^2}+\frac{u^{n + 1}_{j,k + 1}-2u^{n + 1}_{j,k}+u^{n + 1}_{j,k - 1}}{h^2}\right)\\
&=\frac{1}{16h^2}(\delta^2_xu^{n+\frac{1}{2}}_{j,k}+\delta^2_yu^{n + 1}_{j,k})
\end{align*}
\]
用python