# Second derivative using Fornberg finite difference method

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')
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.

• 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 May 20, 2019 at 19:24
• Have you checked this post? scicomp.stackexchange.com/q/11249/9667 May 22, 2019 at 15:58
• 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. May 28, 2019 at 5:16

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