有限元方法中的数值技术：追赶法求解三对角方程
计算原理
A=(b1c1e1b2c2⋱⋱⋱en−2bn−1cn−1en−1bn) \mathbf{A} = \left( \begin{array}{cccccc} b_{1} & c_{1} & & & & \\ e_{1} & b_{2} & c_{2} & & & \\& \ddots & \ddots & \ddots & & \\& & e_{n-2} & b_{n-1} & c_{n-1} & \\& & & e_{n-1} & b_{n} & \\ \end{array} \right) A=b1e1c1b2⋱c2⋱en−2⋱bn−1en−1cn−1bn
对于上述的三对角矩阵，可在Doolittle分解的基础上使用“追赶法”加速矩阵的分解。
L=(1l21l3⋱⋱1ln1),U=(u1d1u2d2⋱⋱un−1dn−1un) \mathbf { L } = \left( \begin{array} { c c c c c c } { 1 } & & & & \\ { l _ { 2 } } & { 1 } & & & \\ & { l _ { 3 } } & { \ddots } & & \\ & & { \ddots } & { 1 } & \\ & & & { l _ { n } } & { 1 } \end{array} \right) , \mathbf { U } = \left( \begin{array} { c c c c c c } { u _ { 1 } } & { d _ { 1 } } & & & \\ & { u _ { 2 } } & { d _ { 2 } } & & \\ & & { \ddots } & { \ddots } & \\ & & & { u _ { n - 1 } } & { d _ { n - 1 } } \\ & & & & { u _ { n } } \end{array} \right) L=1l21l3⋱⋱1ln1,U=u1d1u2d2⋱⋱un−1dn−1un
di=ci,i=1,2,⋯ ,n−1u1=b1
\begin{array} { l } { { d _ { i } = c _ { i } , i = 1 , 2 , \cdots , n - 1 } } \\ { { u _ { 1 } = b _ { 1 } } } \end{array}
di=ci,i=1,2,⋯,n−1u1=b1
li=eiui−1ui=bi−lici−1}i=2,3,⋯ ,n
\left. \begin{array} { l } { \displaystyle l _ { i } = \frac { e _ { i } } { u _ { i - 1 } } } \\ { \displaystyle u _ { i } = b _ { i } - l _ { i } c _ { i - 1 } } \end{array} \right\} i = 2 , 3 , \cdots , n
li=ui−1eiui=bi−lici−1}i=2,3,⋯,n
此时原方程Ax=f\mathbf{A}\mathbf{x}=\mathbf{f}Ax=f变为：
{L y=fU x=y \left\{ \begin{array}{l} \mathbf{L\, y} = \mathbf{f} \\ \mathbf{U\, x} = \mathbf{y} \end{array} \right. {Ly=fUx=y
{y1=f1yi=fi−liyi−1,i=2,3,⋯ ,n \left\{ \begin{array}{l} y_{1} = f_{1} \\ y_{i} = f_{i} - l_{i} y_{i-1}, \quad i = 2,3,\cdots, n \end{array} \right. {y1=f1yi=fi−liyi−1,i=2,3,⋯,n
{xn=ynunxi=(yi−cixi+1)ui,i=n−1,n−2,⋯ ,1 \left\{ \begin{array}{l} x_{n} = \dfrac{y_{n}}{u_{n}} \\ x_{i} = \dfrac{\left(y_{i} - c_{i}x_{i+1}\right)}{u_{i}}, \quad i = n-1, n-2, \cdots, 1 \end{array} \right. ⎩⎨⎧xn=unynxi=ui(yi−cixi+1),i=n−1,n−2,⋯,1
Fortran数值实现
subroutine chase(a,f,x,n)
! chase追赶法求解线性方程组(仅适用于三对角方程组)! Ax=finteger, intent(in) :: nreal*8, intent(in) :: a(n,n), f(n)real*8, intent(out) :: x(n)real*8 :: l(2:n), b(n), u(n), d(1:n-1)real*8 :: c(1:n-1),e(2:n)integer :: ireal*8 ::y(n)!--------------将A矩阵复制给向量e,b,c---------do i =1,nb(i) = a(i,i)end dodo i = 1,n-1c(i) = a(i,i+1)end do do i = 2,ne(i) = a(i,i-1)end do!--------------------------------------------do i = 1,nd(i) = c(i)end dou(1) = b(1)do i = 2,nl(i) = e(i)/u(i-1)u(i) = b(i) - l(i)*c(i-1)end do!----------------回带，求y--------------------y(1) = f(1)do i = 2,ny(i) = f(i) - l(i)*y(i-1)end do!----------------回带，求x--------------------x(n) = y(n)/u(n)do i = n-1,1,-1x(i) = (y(i) - c(i)*x(i+1))/u(i)end do
end subroutine chase
数值案例
(2−100−12−100−12−100−12)x=(1001)
\begin{pmatrix} 2 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 2 \end{pmatrix} \mathbf{x} = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 1 \end{pmatrix}
2−100−12−100−12−100−12x=1001
主程序：
program mainimplicit noneinteger, parameter :: n = 4real*8 :: a(n,n), f(n), x(n)integer :: ia = 0.0d0do i = 1, na(i,i) = 2.0d0 ! 主对角线元素if (i < n) thena(i,i+1) = -1.0d0 ! 上对角线元素a(i+1,i) = -1.0d0 ! 下对角线元素end ifend dof = [1.0d0, 0.0d0, 0.0d0, 1.0d0]call chase(a, f, x, n)write(*,'(/A)') 'Solution vector x:'do i = 1, nwrite(*,'("x(",I1,") = ",F12.6)') i, x(i)end doend program main
输出：
Solution vector x:
x(1) = 1.000000
x(2) = 1.000000
x(3) = 1.000000
x(4) = 1.000000
注意：子程序中默认：
implicit real*8(a-z)
参考文献
- 宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.