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.