Finite difference approximation of Jacobian 1

[dibas]

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:

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

Output(forwad diff.)

 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

 

[mikrom]

You can wruite the forward difference formula as following

       ! Forward difference formula
       df(:,j) = (funcv(x) - funcv(xsav)) / h(j)

Then the central difference formula would be something like this:

subroutine fdjac2(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, xph1, xph2, h, h_double, x2, x1

     n = size(x)

     xsav = x
     h = EPS*abs(xsav)

     where (h <= EPS) h = EPS

     xph2 = xsav + h              ! Trick to reduce finite precision error.
     xph1 = xsav - h
     h_double = xph2 - xph1
     !
     x2 = x
     x1 = x
  do j = 1 , n
       x2(j) = xph2(j)
       x1(j) = xph1(j)
       ! Central difference formula
       df(:,j) = (funcv(x2) - funcv(x1)) / h_double(j)
       x2(j) = xsav(j)
       x1(j) = xsav(j)
  end do

end subroutine fdjac2

I got following results

 Point x at which the Jacobian is evaluated
 --------------------------------
    1.00000000000000000       1.00000000000000000               
 --------------------------------
[COLOR=blue]using Forward difference:

Jacobian matrix, df, at the point x above
-----------------------------------------------

-7.3905340209543704 7.3897950414679103
-29.558441211019140 -22.168646169542033

using Central difference:[/code]
Jacobian matrix, df, at the point x above
-----------------------------------------------

-7.3890562097672312 7.3890561358751841
-29.556224494241849 -22.167168358360133
[/code]

 

[mikrom]

Oops, I got following results:

 Point x at which the Jacobian is evaluated
 --------------------------------
    1.00000000000000000       1.00000000000000000               
 --------------------------------

[COLOR=blue]using Forward difference:[/color]
 Jacobian matrix, df, at the point x above
 -----------------------------------------------

    -7.3905340209543704               7.3897950414679103     
    -29.558441211019140              -22.168646169542033     

[COLOR=blue]using Central difference:[/color]
 Jacobian matrix, df, at the point x above
 -----------------------------------------------

    -7.3890562097672312               7.3890561358751841     
    -29.556224494241849              -22.167168358360133

 

[dibas]

Oops, the central difference estimation is amazing!
Thanks indeed, mikrom!!!