有限元方法中的数值技术：行列式、求逆、矩阵方程

行列式

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

$$\det \left(\mathbf{A}\right)=\det \left(\mathbf{L}\right)\det \left(\mathbf{U}\right)=\det \left(\mathbf{L}\right)=\prod _{i=1}^n{l}_{ii}$$

$$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 mainuse matriximplicit noneinteger, parameter :: n = 5real*8 :: a(n,n), l(n,n), u(n,n)real*8 :: detinteger :: ia = 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.0d0do i = 1, ndet = det * l(i,i)end dowrite(*,'(/A,F12.6)') 'Determinant of A = ', detend 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)$，将上式等价为：

$$\{\begin{array}{c}{\mathbf{AX}}_1={\mathbf{B}}_1\\ {\mathbf{AX}}_2={\mathbf{B}}_2\\ ⋮\\ {\mathbf{AX}}_m={\mathbf{B}}_m\end{array}$$

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

Fortran数值实现

subroutine mat_eq(a,b,x,n,m)
! 矩阵方程求解(既可以计算线性方程组Ax=b，也可以计算矩阵方程AX=B)integer, intent(in) :: n,mreal*8, intent(in) :: a(n,n),b(n,m)real*8, intent(out) :: x(n,m)integer :: i,k,id_maxreal*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) = aab(:,n+1:n+m) = b! 列主元高斯消元do k = 1,n-1elmax = abs(ab(k,k))id_max = kdo i = k+1,nif (abs(ab(i,k)) > elmax) thenelmax = ab(i,k)id_max = iend ifend dovtemp1 = ab(k,:)vtemp2 = ab(id_max,:)ab(k,:) = vtemp2ab(id_max,:) = vtemp1! 进行消元do i = k+1,ntemp = ab(i,k)/ab(k,k)ab(i,:) = ab(i,:) - temp*ab(k,:)end doend doa_up(:,:) = ab(1:n,1:n)! 回带do i = 1,mvtemp = ab(:,n+i)call uptri(a_up,vtemp,xtemp,n)! 将计算结果赋值给xx(:,i) = xtempend doend subroutine mat_eq

数值案例

$$\mathbf{A}=\left(\begin{array}{cccc}1.80& 2.88& 2.05& -0.89\\ 525.00& -295.00& -95.00& -380.00\\ 1.58& -2.69& -2.90& -1.04\\ -1.11& -0.66& -0.59& 0.80\end{array}\right),\mathbf{B}=\left(\begin{array}{cc}9.52& 18.47\\ 2435.00& 225.00\\ 0.77& -13.28\\ -6.22& -6.21\end{array}\right)$$

主程序：

program mainuse matriximplicit noneinteger, parameter :: n = 4, m = 2real*8 :: a(n,n), b(n,m), x(n,m)integer :: i! 使用reshape按行填充矩阵A和Ba = 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, nwrite(*,'(2F12.6)') x(i,:)end doend program main

输出：

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

求逆

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

$$\mathbf{AX}=\mathbf{I}$$

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

Fortran数值实现

subroutine inv_mat(a,inv_a,n)
! 矩阵求逆integer, intent(in) :: nreal*8, intent(in) :: a(n,n)real*8, intent(out) :: inv_a(n,n)real*8 :: e(n,n)integer :: ie = 0.0d0! 初始化单位矩阵edo i = 1, ne(i,i) = 1.0d0end docall mat_eq(a, e, inv_a, n, n)end subroutine inv_mat

数值案例

$$\mathbf{A}=\left(\begin{array}{lll}1& 2& -1\\ 1& 1& 2\\ 3& -1& 1\end{array}\right)$$

输出：

 inv_A:0.176471   -0.058824    0.2941180.294118    0.235294   -0.176471-0.235294    0.411765   -0.058824

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

参考文献

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