Second derivative using Fornberg finite difference method

Question

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
  close(10)


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


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


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


  !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
          enddo
          carray(i,0) = -c1*c5*carray(i-1,0)/c2
      endif

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

  c1 = c2



  enddo

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.

Answer 1

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.

Answer 2

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.