使用泰勒级数展开的 Fortran 矩阵的指数

Question

我有一个 n x n 的矩阵，我想用 Fortran 计算 exponential_of_matrix(matrix_name)。有谁知道用泰勒级数展开计算矩阵的指数吗？

Taylor Series Expansion of e^x = 1 + x + x^2/2! + x^3/3! + x^4/4! + .....

我尝试编写一个矩阵乘法子程序，这在编写矩阵指数子程序时可能会有用。

! mat_mul: matrix_a(n,m) matrix_b(m,l) product(n,l)

subroutine mat_mul(matrix_a,matrix_b, product,n,m,l)
    real*8 matrix_a(n,m), matrix_b(m,l), product(n,l)
    integer i,j,k
    do i=1,n
        do j=1,l
            product(i,j) = 0
            do k=1,m
                product(i,j) = product(i,j) + matrix_a(i,k) * matrix_b(k,j)
            end do
        end do
    end do
end subroutine mat_mul

Answer 1

首先，如果你想乘矩阵，你应该真正使用 matmul¹ 而不是手动操作，所以

do i=1,n
  do j=1,l
    product(i,j) = 0
    do k=1,m
      product(i,j) = product(i,j) + matrix_a(i,k) * matrix_b(k,j)
    end do
  end do
end do

变成

product = matmul(matrix_a, matrix_b)

这段代码更清晰，速度也更快。

第二，你确定要用泰勒级数吗？如果矩阵可对角化，则矩阵求幂通常通过 diagonalisation:

计算

将矩阵对角化以获得其特征向量 v_i 和特征值 e_i。
取特征值的指数e^e_i。
构造具有特征向量 v_i 和特征值 e^e_i 的矩阵。

¹ 如果您需要更高的性能，甚至 BLAS/LAPACK 分布。

Answer 2

我用泰勒级数展开写了这个子程序。而不是 mat_mul 子例程，最好使用 MATMUL(matrix_A,matrix_B).

subroutine exponent_matrix(a,expo_mat,n)
    real a(n,n), prod(n,n), expo_mat(n,n), term(n,n)
    do i=1,n
        do j=1,n
            if(i .eq. j) then
                expo_mat(i,j) =1
                term(i,j) = 1
            else
                expo_mat(i,j) = 0
                term(i,j) = 0
            end if
        end do
    end do
    
    do i=1,5000
        call mat_mul(term,a,prod,n,n,n)
        term = prod/i
        expo_mat = expo_mat + term
    end do
end subroutine exponent_matrix

使用泰勒级数展开的 Fortran 矩阵的指数

exponential of a matrix in Fortran using Taylor Series expansion

fortran