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, 2019 at 19:24
  • $\begingroup$ Have you checked this post? scicomp.stackexchange.com/q/11249/9667 $\endgroup$
    – nicoguaro
    May 22, 2019 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, 2019 at 5:16

2 Answers 2


I agree with what @davidhigh suggested: 300 gridpoints is too much.

In the 1988 paper, Fornberg said that

...the order of accuracy is generally $n-m+1$

for which 300 grid points would lead to 299th order of accuracy if I understood you correctly.


It is to be expected that the weights become large, see Figure 3.2-2 in Fornberg's book A Practical Guide to Pseudospectral Methods. It shows the weights of the first derivative on a uniform grid, but I expect the same trends apply also to second derivatives on non-uniform grids.

Having said that, I doubt that it makes much sense for you to make your stencil as large as the domain. If you consider fewer points in your stencil, for example so you get 4th or 6th order, which is usually good enough, the magnitudes of your weights should be much more reasonable.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.