Hi.
I'm implementing a forward finite difference approximation of the Jacobian
matrix (2x2). I have adapted a subroutine (fdjac) from somewhere.The
issue is that this subroutine requires explicit interfaces because
of the deferred shapes of the involved arrays.
1. May I please have tips on how to eliminate the need for this need
of explicit interfaces. I have attached the code below. I'm sorry,
it looks a unnecessarily lengthy!
2. If anyone has an idea of how to do a CENTRAL finite difference implementation of the Jacobian in f90, I would appreciate their
assistance. Central differences give better approximations that are
closer to the analytical ones. I have included results from the
forward difference coded here, as well as an analytical result from
maple.
Some formulas (multivariate finite diffs. let u = u(x_1, x_2) )
Let eps > 0 (eps small)
Forward diff.: du/dx_1 = ( u(x_1 + eps, x_2) - u(x_1, x_2) ) / eps
du/dx_2 = ( u(x_1, x_2 + eps) - u(x_1, x_2) ) / eps
Central diff.: du/dx_1 = ( u(x_1 + eps, x_2) - u(x_1 - eps, x_2) ) / 2 * eps
du/dx_2 = ( u(x_1, x_2 + eps) - u(x_1, x_2 - eps) ) / 2 * eps
Thank you in advance!
Fortran code:
Output(forwad diff.)
I'm implementing a forward finite difference approximation of the Jacobian
matrix (2x2). I have adapted a subroutine (fdjac) from somewhere.The
issue is that this subroutine requires explicit interfaces because
of the deferred shapes of the involved arrays.
1. May I please have tips on how to eliminate the need for this need
of explicit interfaces. I have attached the code below. I'm sorry,
it looks a unnecessarily lengthy!
2. If anyone has an idea of how to do a CENTRAL finite difference implementation of the Jacobian in f90, I would appreciate their
assistance. Central differences give better approximations that are
closer to the analytical ones. I have included results from the
forward difference coded here, as well as an analytical result from
maple.
Some formulas (multivariate finite diffs. let u = u(x_1, x_2) )
Let eps > 0 (eps small)
Forward diff.: du/dx_1 = ( u(x_1 + eps, x_2) - u(x_1, x_2) ) / eps
du/dx_2 = ( u(x_1, x_2 + eps) - u(x_1, x_2) ) / eps
Central diff.: du/dx_1 = ( u(x_1 + eps, x_2) - u(x_1 - eps, x_2) ) / 2 * eps
du/dx_2 = ( u(x_1, x_2 + eps) - u(x_1, x_2 - eps) ) / 2 * eps
Thank you in advance!
Fortran code:
Code:
!-------------------------------------------------------------------!
! The program main tests the subroutine 'fdjac'.
!
! Description of fdjac and funcv:
!
! Purpose: fdjac computes forward-difference approximation to
! Jacobian.
! Inputs:
! x (length N) - point at which the Jacobian is to be evaluated
!
! fvec (length N) - vector of function values at the point x
!
! funcv(x) - fixed-name, user-supplied routine that returns
! the vector of functions at x
! Output:
! df - N × N Jacobian matrix
!
! Parameter:
! EPS - small perturbation in x.
!
!------------------------------------------------------------------!
program main
implicit none
integer, parameter :: m = 2
double precision, dimension(m) :: u
interface
function funcv(x)
implicit none
double precision, dimension(:), intent(in) :: x
double precision, dimension(size(x)) :: funcv
end function funcv
subroutine fdjac(x,fvec,df)
implicit none
double precision, dimension(:), intent(in) :: fvec
double precision, dimension(:), intent(inout) :: x
double precision, dimension(:,:), intent(out) :: df
double precision :: EPS
integer :: j,n
double precision, dimension(size(x)) :: xsav, xph, h
end subroutine fdjac
end interface
double precision, dimension(m) :: fvec
double precision, dimension(m,m) :: df
u(1) = 1.0
u(2) = 1.0
write(*,*)
write(*,*) "Point x at which the Jacobian is evaluated"
write(*,*) "--------------------------------"
write(*,*) " ", u, " "
write(*,*) "--------------------------------"
! Now call the function 'funcv' to calculate it's value at
! x and store it in (function) 'fvec'. This vector of function
! values is needed when calling the subroutine fdjac
fvec = funcv(u)
call fdjac(u,fvec,df)
write(*,*)
write(*,*) "Jacobian matrix, df, at the point x above"
write(*,*) "-----------------------------------------------"
write(*,*)
write(*,*) " ", df(1,1), " ", df(1,2)
write(*,*) " ", df(2,1), " ", df(2,2)
write(*,*)
end program main
function funcv(x)
implicit none
double precision, dimension(:), intent(in) :: x
double precision, dimension(size(x)) :: funcv
! Want Jacobian matric for the following functions:
funcv(1) = exp( x(1) + x(2) ) * x(2)*x(1) - &
exp( x(1) + x(2) ) * x(1)**2
funcv(2) = 1 - exp( x(1) + x(2) ) * x(2) * x(1) - &
exp( x(1) + x(2) ) * x(1)
end function funcv
subroutine fdjac(x,fvec,df)
implicit none
double precision, dimension(:), intent(in) :: fvec
double precision, dimension(:), intent(inout) :: x
double precision, dimension(:,:), intent(out) :: df
interface
function funcv(x)
implicit none
double precision, dimension(:), intent(in) :: x
double precision, dimension(size(x)) :: funcv
end function funcv
end interface
double precision, parameter :: EPS = 1.0e-04
integer :: j,n
double precision, dimension(size(x)) :: xsav, xph, h
n = size(x)
xsav = x
h = EPS*abs(xsav)
where (h <= EPS) h = EPS
xph = xsav + h ! Trick to reduce finite precision error.
h = xph - xsav
do j = 1 , n
x(j) = xph(j)
df(:,j) = (funcv(x) - fvec(:)) / h(j) ! Forward difference formula
x(j) = xsav(j)
end do
end subroutine fdjacOutput(forwad diff.)
Code:
Point x at which the Jacobian is evaluated
--------------------------------
1.0000000000000000 1.0000000000000000
--------------------------------
Jacobian matrix, df, at the point x above
-----------------------------------------------
-7.3905340209540542 7.3897950414682274
-29.558441211019140 -22.168646169559796
maple output:
−7.389056101 7.389056099
−29.55622440 −22.16716830