追赶法求解三对角矩阵

理论推导

① 用$a_{ij}$表示

从如下的出发

$$\left[\begin{array}{ccccc}{a}_{12}& a_{13}& 0& \cdots & 0\\ a_{21}& a_{22}& a_{23}& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& \cdots & a_{N-1,1}& a_{N-1,2}& a_{N-1,3}\\ 0& \cdots & 0& a_{N,1}& a_{N,2}\end{array}\right]\left[\begin{array}{c}{u}_1^{k+1}\\ u_2^{k+1}\\ \cdots \\ u_{N-1}^{k+1}\\ u_N^{k+1}\end{array}\right]=\left[\begin{array}{c}{R}_1\\ R_2\\ \cdots \\ R_{N-1}\\ R_N\end{array}\right]$$

可以进一步整理成

$$\begin{array}{rl}& \left[\begin{array}{ccccc}{a}_{12}& a_{13}& 0& \cdots & 0\\ 0& a_{22}^{\prime }& a_{23}& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& \cdots & 0& a_{N-1,2}^{\prime }& a_{N-1,3}\\ 0& \cdots & 0& 0& a_{N,2}^{\prime }\end{array}\right]\left[\begin{array}{c}{u}_1^{k+1}\\ u_2^{k+1}\\ \cdots \\ u_{N-1}^{k+1}\\ u_N^{k+1}\end{array}\right]=\left[\begin{array}{c}{R}_1^{\prime }\\ R_2^{\prime }\\ \cdots \\ R_{N-1}^{\prime }\\ R_N^{\prime }\end{array}\right]\end{array}$$

其中$R_1^{\prime }=R_1$

$$\{\begin{array}{l}{a}_{i2}^{\prime }=a_{i2}-\frac{a_{i1}{a}_{i-1,3}}{a_{i-1,2}^{\prime }}\\ R_i^{\prime }=R_i-\frac{a_{i1}{R}_{i-1}^{\prime }}{a_{i-1,2}^{\prime }}\end{array}(2\le i\le N,a_{12}^{\prime }=a_{12},a_{13}^{\prime }=a_{13})$$

基于整理后的矩阵，可以得到

$$\Rightarrow \{\begin{array}{l}{u}_N^{k+1}=\frac{R_N^{\prime }}{a_{N,2}^{\prime }}\\ u_i^{k+1}=\frac{R_i^{\prime }-a_{i3}{u}_{i+1}^{k+1}}{a_{i2}^{\prime }}(1\le i\le N-1)\end{array}$$

也就是说，先通过最后一行计算得到$u_N^{k+1}$，然后再从最后一行往上计算。

② 用$a_i,b_i,c_i$表示

从如下的三对角矩阵出发，$AU=R$

$$\left[\begin{array}{ccccc}{b}_1& c_1& 0& \cdots & 0\\ a_2& b_2& c_2& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& \cdots & a_{N-1}& b_{N-1}& c_{N-1}\\ 0& \cdots & 0& a_N& b_N\end{array}\right]\left[\begin{array}{c}{u}_1^{k+1}\\ u_2^{k+1}\\ \cdots \\ u_{N-1}^{k+1}\\ u_N^{k+1}\end{array}\right]=\left[\begin{array}{c}{R}_1\\ R_2\\ \cdots \\ R_{N-1}\\ R_N\end{array}\right]$$

其中$b_i$对角元，$a_i$左旁元，$c_i$右旁元. 可以进一步整理成

$$\begin{array}{rl}& \left[\begin{array}{ccccc}{b}_1& c_1& 0& \cdots & 0\\ 0& b_2^{\prime }& c_2& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& \cdots & 0& b_{N-1}^{\prime }& c_{N-1}\\ 0& \cdots & 0& 0& b_N^{\prime }\end{array}\right]\left[\begin{array}{c}{u}_1^{k+1}\\ u_2^{k+1}\\ \cdots \\ u_{N-1}^{k+1}\\ u_N^{k+1}\end{array}\right]=\left[\begin{array}{c}{R}_1^{\prime }\\ R_2^{\prime }\\ \cdots \\ R_{N-1}^{\prime }\\ R_N^{\prime }\end{array}\right]\end{array}$$

其中$R_1^{\prime }=R_1$，且有

$$\{\begin{array}{l}{b}_i^{\prime }=b_i-\frac{a_i{c}_{i-1}}{b_{i-1}^{\prime }}\\ R_i^{\prime }=R_i-\frac{a_i{R}_{i-1}^{\prime }}{b_{i-1}^{\prime }}\end{array}(2\le i\le N,b_1^{\prime }=b_1,c_1^{\prime }=c_1)$$

基于整理后的矩阵，可以得到

$$\Rightarrow \{\begin{array}{l}{u}_N^{k+1}=\frac{R_N^{\prime }}{b_N^{\prime }}\\ u_i^{k+1}=\frac{R_i^{\prime }-c_i{u}_{i+1}^{k+1}}{b_i^{\prime }}(1\le i\le N-1)\end{array}$$

也就是说，先通过最后一行计算得到$u_N^{k+1}$，然后再从最后一行往上计算.

伪程序/流程

求解$AU=R$

$a$左旁元 → a(ny)

$b$对角元 → b(ny)

$c$右旁元 → c(ny)

$R$ → f(ny)

$u_i^{k+1}$ → u(i,n+1)，维度是u(ny,nt)

ny：实空间网格数， nt：时间空间网格数

• 计算$b_i^{\prime }=b_i-\frac{a_i{c}_{i-1}}{b_{i-1}^{\prime }}$

        M(1)=b(1)
        do i=2,ny-1
          L(i)=a(i)/M(i-1)
        M(i)=b(i)-L(i)*c(i-1)
        end do

其中M(i)即为$b_i^{\prime }$

• 计算$R$

        do i = 1,ny
         f(i)=R(i)
        end do

一般来说，R(i)与矩阵元素$a,b,c$以及$u$均有关系，这里为了方便理解，以伪代码的形式展示，故写成以上形式，具体实践见CN差分格式求解偏微分方程

• 计算$R_i^{\prime }=R_i-\frac{a_i{R}_{i-1}^{\prime }}{b_{i-1}^{\prime }}$

        y(1)=f(1)
        do i=2,ny
         y(i)=f(i)-L(i)*y(i-1)
        end do

注意， L(i)=a(i)/M(i-1)

• 计算$u_N^{k+1}=\frac{R_N^{\prime }}{b_N^{\prime }},u_i^{k+1}=\frac{R_i^{\prime }-c_i{u}_{i+1}^{k+1}}{b_i^{\prime }}(1\le i\le N-1)$

        u(ny,n+1)=y(ny)/M(ny)
        do i=ny-1,1,-1
          u(i,n+1)=(y(i)-c(i)*u(i+1,n+1))/M(i)
        end do

至此可以计算出n+1时刻的整个空间网格上的$u$，在时间循环里不断调用，就可以计算出整个时间网格上的$u$