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