Search code examples
fortran

exponential of a matrix in Fortran using Taylor Series expansion


I am having a matrix of n x n and I want to calculate exponential_of_matrix(matrix_name) in Fortran. Is there anyone who knows to calculate the exponential of a matrix using Taylor Series Expansion?

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

I tried to write a matrix multiplication subroutine which may be useful during writing matrix exponential subroutine.

! 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

Solution

  • I have written this subroutine using Taylor Series Expansion. Instead of mat_mul subroutine, its better to use 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