# 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 – Ian Bush May 20 at 19:24
• Have you checked this post? scicomp.stackexchange.com/q/11249/9667 – nicoguaro May 22 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. – davidhigh May 28 at 5:16