I have some discrete data, non-equispaced in x, y=f(x).

I want to use a numerical finite difference method to calculate the second derivatives of y, at some point.

I am using the Fornberg method, which is well described here and here, and working in Fortran.

Program fornberg

  implicit none
  integer, parameter :: dp = selected_real_kind(32,307)
  integer, parameter :: Nrows =1d4 
  integer (kind=dp) :: j, counts,m,k, N, i, nn, mn
  real(kind=dp), dimension(Nrows,2) :: SomeData
  real(kind=dp), dimension(:,:), allocatable ::carray
  real(kind=dp), dimension(:), allocatable :: xdata,ydata
  real(kind=dp) :: c1,c2,c3,c4,c5, z

  !Load the data from an unformatted binary file

  open(unit = 10, file='example.dat', form = 'unformatted')
  read(10) SomeData

  !Count the number of non-zero rows
  do j=1,Nrows
  if (SomeData(j,1) .EQ. 0) then
  counts = j-1

  !Allocate the x and y arrays to hold the data with indexing starting at 0
  ALLOCATE(xdata(0:counts-1)) !specific zero indexing

  !Populate the arrays
  do j = 1,counts
  xdata(j-1) = SomeData(j,1)
  ydata(j-1)= SomeData(j,2)

  !Define the point at which we want the derivative evaluated
  z = xdata(0)

  !Length of data array
  N = counts
  nn = N - 1

  !Define the maximum order of the derivative
  m = 2

  !Set up zeroes array
  ALLOCATE(carray(0:N-1, 0:m)) !zero indexing

  !Determine the weights via the Fornberg algorithm

  c1 = 1.0_dp
  c4 = xdata(0) - z
  carray = 0.0_dp
  carray(0,0) = 1.0_dp

  do i = 1,nn

  mn = min(i,m)

  c2 = 1.0_dp
  c5 = c4
  c4 = xdata(i) - z

  do j = 0, i-1

      c3 = xdata(i) - xdata(j)
      c2 = c2*c3

      if ( j .EQ. i-1) then

          do k=mn,1,-1
              carray(i,k) = c1*(k*carray(i-1,k-1) - c5*carray(i-1,k))/c2
          carray(i,0) = -c1*c5*carray(i-1,0)/c2

      do k = mn,1,-1
          carray(j,k) = (c4*carray(j,k) - k*carray(j,k-1))/c3
      carray(j,0) = c4*carray(j,0)/c3

  c1 = c2


end program fornberg

My problem is that the weights for the second derivatives quickly become huge.

Looking at the code, this is a direct consequence of the c2 = c2*c3 command. For c3 > 1 and a large number of iterations (the dataset is ~300 rows), I am confused about how the weights could ever be 'reasonable'.

Any guidance would be greatly appreciated. I can also provide the dataset if necessary.

  • 1
    $\begingroup$ What you want is c1/c2 - a quick glance suggests that while both those numbers will get huge, the ratio will stay sensible. Try rewriting the code to work in terms of that ratio $\endgroup$ – Ian Bush May 20 at 19:24
  • $\begingroup$ Have you checked this post? scicomp.stackexchange.com/q/11249/9667 $\endgroup$ – nicoguaro May 22 at 15:58
  • $\begingroup$ Are you working on an equally spaced grid? If so, the problem is exponentially ill-conditioned and 300 gridpoints is much too much. Better use Chebyshev points then. $\endgroup$ – davidhigh May 28 at 5:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.