行列式、求逆、矩阵方程¶

约 525 个字 130 行代码预计阅读时间 3 分钟

行列式的计算¶

在矩阵完成$LU$分解后，可以直接利用单位三角矩阵进行计算，以“Crout”分解为例 $$ \mathbf { A } = \mathbf { L U } $$

\[ \operatorname* { d e t } \left( \mathbf { A } \right) = \operatorname* { d e t } \left( \mathbf { L } \right) \operatorname* { d e t } \left( \mathbf { U } \right) = \operatorname* { d e t } \left( \mathbf { L } \right) = \prod _ { i = 1 } ^ { n } l _ { i i } \]

\[ A = \left( \begin{array}{ccccc} 2.7486 & 1.9022 & 3.8958 & 0.0595 & 2.6427 \\ 4.5860 & 2.8391 & 4.6701 & 1.6856 & 0.8282 \\ 1.4292 & 0.3793 & 0.6495 & 0.8109 & 3.0099 \\ 3.7860 & 0.2698 & 2.8441 & 3.9714 & 1.3149 \\ 3.7686 & 2.6540 & 2.3470 & 1.5561 & 3.2704 \end{array} \right) \]

主程序：

program main
    use matrix
    implicit none

    integer, parameter :: n = 5
    real*8 :: a(n,n), l(n,n), u(n,n)
    real*8 :: det
    integer :: i

    a = reshape([2.7486d0, 1.9022d0, 3.8958d0, 0.0595d0, 2.6427d0, &
                 4.5860d0, 2.8391d0, 4.6701d0, 1.6856d0, 0.8282d0, &
                 1.4292d0, 0.3793d0, 0.6495d0, 0.8109d0, 3.0099d0, &
                 3.7860d0, 0.2698d0, 2.8441d0, 3.9714d0, 1.3149d0, &
                 3.7686d0, 2.6540d0, 2.3470d0, 1.5561d0, 3.2704d0], &
                [n, n], order=[2,1])

    call crout(a, l, u, n)

    det = 1.0d0
    do i = 1, n
        det = det * l(i,i)
    end do

    write(*,'(/A,F12.6)') 'Determinant of A = ', det

end program main

输出：

Determinant of A =   -22.887183

矩阵方程¶

$$ \mathbf { A } _ { [ N , N ] } \mathbf { X } _ { [ N , M ] } = \mathbf { B } _ { [ N , M ] } $$ 记$\mathbf { X }$的第$i$列向量为$\mathbf { X }_{i}(i=1,m)$，$\mathbf { B }$的第$i$列向量为$\mathbf { B }_{i}(i=1,m)$，将上式等价为：

\[ \left\{ \begin{array} { c } { \mathbf { A X } _ { 1 } = \mathbf { B } _ { 1 } } \\ { \mathbf { A X } _ { 2 } = \mathbf { B } _ { 2 } } \\ { \vdots } \\ { \mathbf { A X } _ { m } = \mathbf { B } _ { m } } \end{array} \right. \]

但在实际计算时不应该分别求解，如果是这样的话就造成计算机资源极大浪费，而应该是对所有向量一次选主元，然后分别回代。

Fortran数值实现¶

subroutine mat_eq(a,b,x,n,m)
! 矩阵方程求解(既可以计算线性方程组Ax=b，也可以计算矩阵方程AX=B)
    integer, intent(in) :: n,m
    real*8, intent(in) :: a(n,n),b(n,m)
    real*8, intent(out) :: x(n,m)

    integer :: i,k,id_max
    real*8 :: a_up(n,n),b_up(n,m)
    real*8 :: ab(n,n+m)
    real*8 :: vtemp1(n+m),vtemp2(n+m)
    real*8 :: vtemp(n),xtemp(n)

    ab(1:n,1:n) = a
    ab(:,n+1:n+m) = b

    ! 列主元高斯消元
    do k = 1,n-1
        elmax = abs(ab(k,k))
        id_max = k
        do i = k+1,n
            if (abs(ab(i,k)) > elmax) then
                elmax = ab(i,k)
                id_max = i
            end if
        end do
        vtemp1 = ab(k,:)
        vtemp2 = ab(id_max,:)

        ab(k,:) = vtemp2
        ab(id_max,:) = vtemp1
        ! 进行消元
        do i = k+1,n
            temp = ab(i,k)/ab(k,k)
            ab(i,:) = ab(i,:) - temp*ab(k,:)
        end do
    end do

    a_up(:,:) = ab(1:n,1:n)
    ! 回带
    do i = 1,m
        vtemp = ab(:,n+i)
        call uptri(a_up,vtemp,xtemp,n)
        ! 将计算结果赋值给x
        x(:,i) = xtemp
    end do

end subroutine mat_eq

数值案例¶

\[ \mathbf { A } = \left( \begin{array} { c c c c } { 1 . 8 0 } & { 2 . 8 8 } & { 2 . 0 5 } & { - 0 . 8 9 } \\ { 5 2 5 . 0 0 } & { - 2 9 5 . 0 0 } & { - 9 5 . 0 0 } & { - 3 8 0 . 0 0 } \\ { 1 . 5 8 } & { - 2 . 6 9 } & { - 2 . 9 0 } & { - 1 . 0 4 } \\ { - 1 . 1 1 } & { - 0 . 6 6 } & { - 0 . 5 9 } & { 0 . 8 0 } \end{array} \right) , \ \ \mathbf { B } = \left( \begin{array} { c c c } { 9 . 5 2 } & { 1 8 . 4 7 } \\ { 2 4 3 5 . 0 0 } & { 2 2 5 . 0 0 } \\ { 0 . 7 7 } & { - 1 3 . 2 8 } \\ { - 6 . 2 2 } & { - 6 . 2 1 } \end{array} \right) \]

主程序：

program main
    use matrix
    implicit none

    integer, parameter :: n = 4, m = 2
    real*8 :: a(n,n), b(n,m), x(n,m)
    integer :: i

    ! 使用reshape按行填充矩阵A和B
    a = reshape([1.80d0, 2.88d0, 2.05d0, -0.89d0, &
                 525.00d0, -295.00d0, -95.00d0, -380.00d0, &
                 1.58d0, -2.69d0, -2.90d0, -1.04d0, &
                 -1.11d0, -0.66d0, -0.59d0, 0.80d0], &
                [n, n], order=[2,1])

    b = reshape([9.52d0, 18.47d0, &
                 2435.00d0, 225.00d0, &
                 0.77d0, -13.28d0, &
                 -6.22d0, -6.21d0], &
                [n, m], order=[2,1])

    call mat_eq(a, b, x, n, m)

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

end program main

输出：

Solution X:
    1.000000    3.000000
   -1.000000    2.000000
    3.000000    4.000000
   -5.000000    1.000000

求逆¶

设矩阵$\mathbf { A } _ { \left[ N , N \right] }$为 非奇异矩阵，保证逆矩阵存在，记$\mathbf { A } ^ { - 1 } = \mathbf { X }$，则

\[ \mathbf { A X } \! = \! \mathbf { I } \]

将计算矩阵 A 的逆矩阵转化为计算矩阵方程问题，使用上一节的方法求解即可。

Fortran数值实现¶

subroutine inv_mat(a,inv_a,n)
! 矩阵求逆
    integer, intent(in) :: n
    real*8, intent(in) :: a(n,n)
    real*8, intent(out) :: inv_a(n,n)

    real*8 :: e(n,n)
    integer :: i

    e = 0.0d0
    ! 初始化单位矩阵e
    do i = 1, n
        e(i,i) = 1.0d0
    end do

    call mat_eq(a, e, inv_a, n, n)

end subroutine inv_mat

数值案例¶

\[ \mathbf { A } = { \left( \begin{array} { l l l } { 1 } & { 2 } & { - 1 } \\ { 1 } & { 1 } & { 2 } \\ { 3 } & { -1 } & { 1 } \end{array} \right) } \]

输出：

 inv_A:
    0.176471   -0.058824    0.294118
    0.294118    0.235294   -0.176471
   -0.235294    0.411765   -0.058824

注意：子程序中默认：implicit real*8(a-z)

参考文献¶

宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.