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

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

## 2 Answers

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.