计算矩阵函数 exp(A)

[地球屋里老师]

module m
contains
!****************************************************************************
! 计算矩阵函数 exp(A)
! author: LI xingwang, Chang''an university, China.
! email : [e-mail]
!****************************************************************************
subroutine matrixFunExp(a, res, eps)
implicit none
real, intent(in):: a(:,:)
real,intent(inout):: res(:,:)
real,optional, intent(in):: eps
real, allocatable:: b(:,:)
real  e, c
integer k
e = merge(eps,1.0e-7,present(eps))
allocate(b,source=a)
res = a + reshape([1,0,0,0,1,0,0,0,1],shape(a))
k = 1; c = 1.0
do
  k = k + 1
  c = c / k
  b = matmul(b,a)
  res = res + c*b
  if(all(abs(c*b)<e)) exit
end do
end subroutine
end module

program test
use m
implicit none
real a(3,3), b(3,3), t
integer i
t = 0.5
a = reshape([-1,1,-4,0,2,0,1,0,3],shape(a))
call matrixFunExp(a*t, b) !计算 exp(A*t)
print*,'数值解'
do i = 1, size(b,1)
  print*,b(i,:)
end do
print*,'解析解'
a(:,1) =[exp(t)-2*t*exp(t),-exp(2*t)+exp(t)+2*t*exp(t),-4*t*exp(t)]
a(:,2) =[0.0,exp(2*t),0.0]
a(:,3) =[t*exp(t),exp(2*t)-exp(t)-t*exp(t),2*t*exp(t)+exp(t)]
do i = 1, size(b,1)
  print*,a(i,:)
end do
print*,'误差'
do i = 1, size(b,1)
  print*,a(i,:)-b(i,:)
end do
pause
end