Jacobian matrix: discrepancy to Matlab

1

Hello community,
while playing around with numerical calculation of Jacobian matrices by the method of finite differences (e.g. central differences), I noticed, that the same nonlinear equation produces a different Jacobian matrix when calculated with Matlab, compared to the one by Fortran, by ca.10%. The method involves estimating the partial derivatives of a function with respect to the variables and uses very small numbers (1e-07) as the increment dx. Basically:

F’=dF/dx= [ F(x+dx)-F(x-dx) ] / 2*dx
Since there are very small numbers involved in the evaluation of the nominator and also used in the denominator, could it be that Fortran comes to different results? I used double precision numbers in both cases.

Thanks

2

Hi @lanast, thanks for the question. Can you please post your Fortran and Matlab codes? Let’s have a look.

1 Like
3

@Ianast,

Can you please confirm you have checked your Fortran code, especially your function evaluations, using suitable and small unit tests and ensured the oft overlooked items such as precision and arithmetic of the floating-point calculations are all consistent with your needs?

You can review the trivial example below to see the different possibilities.

   blk1: block
      real :: x, dx, dfdx
      print *, "Block 1: use default precision and mixed-mode arithmetic"
      x = 0.4
      dx = sqrt( epsilon(dx) )
      dfdx = ( sin(x+dx) - sin(x-dx) )/ (2*dx)
      print *, "f:= sin(x); at x=0.4, dfdx =  ", dfdx
      print *, "analytical solution: ", cos(x)
      print "(g0,1x,1pg10.2)", "% error", (dfdx/cos(x) - 1.0)*100.0
   end block blk1
   print *
   blk2: block
      integer, parameter :: wp = selected_real_kind( p=12 )
      real(kind=wp) :: x, dx, dfdx
      print *, "Block 2: use defined precision of ", precision(x), " and consistent arithmetic"
      x = 0.4_wp
      dx = sqrt( epsilon(dx) )
      dfdx = ( sin(x+dx) - sin(x-dx) )/ (2.0_wp*dx)
      print *, "f:= sin(x); at x=0.4, dfdx =  ", dfdx
      print *, "analytical solution: ", cos(x)
      print "(g0,1x,1pg10.2)", "% error", (dfdx/cos(x) - 1.0_wp)*100.0_wp
   end block blk2
end

Upon execution of the program, the expected output is

Block 1: use default precision and mixed-mode arithmetic
f:= sin(x); at x=0.4, dfdx = 0.921042860
analytical solution: 0.921060979
% error -1.97E-03

Block 2: use defined precision of 15 and consistent arithmetic
f:= sin(x); at x=0.4, dfdx = 0.92106099426746368
analytical solution: 0.92106099400288510
% error 2.87E-08

As you know, such small unit tests can help you considerably with debugging the differences.

3 Likes
4

Hi,

I’ve used finite different schemes (central differences or higher order form) with fortran and I don’t have inconsistencies with respect to analytical values (when available).
Unless, you have some inconsistencies with the real kind in your code, you can try with a larger dx (10^-4). With this value and for x value around one, the error (fana’(x)-fnum’(x)) is not larger then 10^-7 for all intrinsic functions.

2 Likes
5

Thank you all for the suggestions. I will have to do some research wheter the epsilon is suitable here. The code is rather clumsy to be posted here, since subroutines are called within. I noticed that by varying the dx increment of the finite differences, either the matrix columns corresponding to displacement or the columns corresponding to velocities agree with Matlab, but not at the same time simultaneously. So with dx=1e-7 the odd columns are identical and the even not, while with dx=1e-4 it is vice ersa.I guess it has something to do with the scaling of the variables in my simulation…