追赶法求解三对角方程¶

约 358 个字 73 行代码预计阅读时间 2 分钟

计算原理¶

\[ \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) \]

对于上述的三对角矩阵，可在Doolittle分解的基础上使用“追赶法”加速矩阵的分解。

\[ \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) \]

\[ \begin{array} { l } { { d _ { i } = c _ { i } , i = 1 , 2 , \cdots , n - 1 } } \\ { { u _ { 1 } = b _ { 1 } } } \end{array} $$ $$ \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 \]

此时原方程$\mathbf{A}\mathbf{x}=\mathbf{f}$变为：

\[ \left\{ \begin{array}{l} \mathbf{L\, y} = \mathbf{f} \\ \mathbf{U\, x} = \mathbf{y} \end{array} \right. \]

\[ \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. \]

\[ \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. \]

Fortran数值实现¶

subroutine chase(a,f,x,n)
! chase追赶法求解线性方程组(仅适用于三对角方程组)
    ! Ax=f
    integer, intent(in) :: n
    real*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 :: i
    real*8 ::y(n)

    !--------------将A矩阵复制给向量e,b,c---------
    do i =1,n
        b(i) = a(i,i)
    end do
    do i = 1,n-1
        c(i) = a(i,i+1)
    end do 
    do i = 2,n
        e(i) = a(i,i-1)
    end do
    !--------------------------------------------
    do i = 1,n
        d(i) = c(i)
    end do
    u(1) = b(1)
    do i = 2,n
        l(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,n
        y(i) = f(i) - l(i)*y(i-1)
    end do
    !----------------回带，求x--------------------
    x(n) = y(n)/u(n)
    do i = n-1,1,-1
        x(i) = (y(i) - c(i)*x(i+1))/u(i)
    end do
end subroutine chase

数值案例¶

$$ \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} $$ 主程序：

program main
    implicit none

    integer, parameter :: n = 4
    real*8 :: a(n,n), f(n), x(n)
    integer :: i

    a = 0.0d0
    do i = 1, n
        a(i,i) = 2.0d0  ! 主对角线元素
        if (i < n) then
            a(i,i+1) = -1.0d0  ! 上对角线元素
            a(i+1,i) = -1.0d0  ! 下对角线元素
        end if
    end do

    f = [1.0d0, 0.0d0, 0.0d0, 1.0d0]

    call chase(a, f, x, n)

    write(*,'(/A)') 'Solution vector x:'
    do i = 1, n
        write(*,'("x(",I1,") = ",F12.6)') i, x(i)
    end do

end 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.