First Derivative

[Nazmul_I_1]

Hi,

I have some data (z,h) which reprents sinusoidal curve. Now, I would like to calculate first derivative dh/dz at each point (z,h).

I am using Intel Fortran Compiler with IMSL.

Can you please help me in this regard?

Cheers

[mecej4]

Do your data satisfy the Nyquist sampling condition? If not, the task may be impossible. Do you know the frequency of the sinusoidal function? How many data points do you have, and are they subject to errors?

[Anthony_Richards]

Use a cubic spline routine to fit to your data points and get a set of cubic spline interpolation coefficients.

Then compute the function values at two points very close to  but on either side of the value where you want the differential coefficient and compute an approximation to it by taking the difference between the two function values and dividing by the difference between the two abscissa values.

Better is to differentiate directly the fitted function polynomial function, if possible, and then evaluate the DC at the abscissa value you want.

Generally, a piecewise cubic function Y(X) = A0+A1*X+A2*X^2+A3*X^3  valid over  a range X1<X<X2 will differentiate to

DYDX(X)=A1+2A2*X+3*A3*X^2

So when you know the A-coefficients of the function fitted over the range X1<X<X2, you can find the DC within the range by using them too.

[jimdempseyatthecove]

Nazmul,

>>I have some data (z,h)...

How is this stored?

From your description, one might assume you have a series of N sample points, each sample point contains a Time value (z) and a measured value (h). As mecej4 asks, does the data contain errors, and if so, what is the nature of the errors. Based on this assumption an no error in data.

integer, parameter :: z=1
integer, parameter :: h=2
integer :: nPoints
real, allocatable :: data(:, :) ! 2, nPoints
...
nPoints = get_nPoints()
allocate data(2,nPoints)
call get_data(data, nPoints)
! if necessary
call sort_data(data, nPoints, z) ! assure ascending z
...

real function get_dhdz(data, nPoints, at_some_z)
  implicit none
  integer :: nPoints
  real :: data(2,nPoints)
  real :: at_some_z
  integer :: i
  if(nPoints .lt. 2) stop "two few of points"
  do i=1,nPoints-1
    if(data(z,i) .eq. at_some_z) then
      if(i .eq. 1) then
        get_dhdz = (data(h,i+1)-data(h,i)) / (data(z,i+1)-data(z,i))
      else
        get_dhdz = (data(h,i+1)-data(h,i-1)) / (data(z,i+1)-data(z,i-1))
      endif
      return
    endif
    if(data(z,i) .gt. at_some_z) then
      if(i .eq. 1) then
        get_dhdz = (data(h,i+1)-data(h,i)) / (data(z,i+1)-data(z,i))
      else
        get_dhdz = (data(h,i)-data(h,i-1)) / (data(z,i)-data(z,i-1))
      endif
      return
    endif
  end do
  i = nPoints
  get_dhdz = (data(h,i)-data(h,i-1)) / (data(z,i)-data(z,i-1))
end function get_dhdz

If the data has errors, then you may need to perform a curve fit and/or other filtering techniques on the data. It might be easiest to conceive of producing data_filtered from data, then apply the simple routine above. Note, data_filtered could be designed to be much larger than data, containing interpolated points that fit the curve that meets your requirements.

Jim Dempsey

[TimP]

The suggestions by Anthony and Jim have much different behavior in presence of noisy data.  Chebyshev fits may be useful but seem beyond scope of this forum.

[Nazmul_I_]

Thanks all for your valuable comments.

Regarding mecej4 comment: I am not sure about ''Nyquist sampling condition''. But I know the frequency of the sinusoidal function. Also, I have about 100 data points (z, h(z)) with no error in the data.

Thanks Anthony for giving your nice idea.

Thanks Jim for proving the codes.

Thanks Tim for your comments.

[jimdempseyatthecove]

If you do not have sufficient number of data points... per half wave length then you really should take Anthony's suggestion and determine a suitable curve that fits the points and then use the curve for interpolation.

An alternative might be to use a modified version of what I listed above and for "at_some_z", determine the slope at the point of data just below at_some_z and just after at_some_z. Then use the distance between the two measured points where at_some_z lies, and use this to interpolate. You will have to do a little more work when the signs of the slopes of the bounding points are different. This should be QED.

Jim Dempsey

[Nazmul_I_]

Thanks Jim for your suggestion.

[John_Campbell]

Nazmul,

If you know the frequency, it would be interesting to calculate the error in your sample points.

The issue of noisy data and how to manage it is a very interesting topic and with the benefit of modern computers it can be interesting to review the different numerical approaches that were developed before computers were available.

The use of Chebyshev or Bessel interpolation formula can be very useful. Higher order functions have their place when the sample points are accurate, but all real numbers have some error and any error can result in unusual results.

Where noise is expected, I find a least squares fit is a useful approach as it also provides some error analysis.

Another approach I have successfully used is a 4 point Lagrange interpolation formula, based on the shape function approach used in finite element formulations. If the sample points are evenly spaced, the shape functions can be easily defined and the derivative is easy to derive. This is much easier than the use of other high order interpolation approaches.

[Anthony_Richards]

For a problem as straightforward as this, my first call would actually be to utilise my MathCad 15 package.
The fitting of a sinusoudal function Y(X)=A+B*sin(C*X+D) to determine best-fit  A,B,C and D values can be done in two lines of code.
Once the coefficients are available, the differential coefficient function Y(x)=C*B*cos(C*X+D) can then immediately be used for any value of X. MathCad's cubic spline interpolation capabilities are also very easy and compact to use for simple problems that do not need you to go through the whole Fortran program + data input route. I use MathCad regularly to attack an algorithm and only switch to Fortran if intensive computation is required.

[Nazmul_I_]

Thanks John for your different suggestions.

Thanks Anthony for your suggestion regarding MathCad. Actually, I have four different curves - sinusoidal, bell shaped, triangular and rectangular. Am I able to evaluate first derivative in all cases by MatCad?

[Anthony_Richards]

Forget about MathCad. You are going to have to do some work.
The task is easier if you have equi-spaced data points.

The simplest way of obtaining the first derivative at a point x=x0+p*h (0< or = p< or = 1) from a data set consisting of
points uniformly spaced at intervals h in x is to first appreciate how to interpolate a value for Y=f(x) at x=x0+p*h from
three values of Y =f(x) at x=x0-h, x=x0 and x=x0+h. This method uses the simple Lagrange 3-point interpolation formula, namely

f(x0+p*h)=0.5*p*(p-1)f(x0-h) + (1-p^2)*f(x0) + 0.5*p*(p+1)*f(x0+h)

Differentiating this wrt p gives you a formula for the first differential coefficient, namely:

df(x0*p*h)/dp  = 0.5*(2*p-1)f(x0-h) - 2*p*f(x0) + 0.5*(2*p+1)*f(x0+h)

So, you take your X-value, determine which, in your uniform list of x-values, is the nearest x-value to it (x0) and the
two values on either side of it, (x0-h) and (x0+h) then use the corresponding Y-values at the three points x0-h, x0
and x0+h, together with the quantity p= (X-x0)/h and use them in the formula

dY/dx= 0.5*(2*p-1)f(x0-h) - 2*p*f(x0) + 0.5*(2*p+1)*f(x0+h)

It will be a nice programming exercise for you.

Questions about accuracy arise wherever the interpolation formula is likely to be inaccurate because the
function may be changing more rapidly than the density of points can describe or follow, but if you have well-behaved
'smooth' functions with a reasonable uniform density of x-values, you should be OK.

Let us know how you get on!

[Nazmul_I_]

Hi Anthony,

Many thanks for your advice. However, the data are not uniformly spaced and I have approximately 100 data points. But I have tested the calculation of first derivative with non-uniform Lagrange 3 point interpolating formula which are merely close to the exact value. Also, as you and Jim suggested, I have calculated the first derivative by using IMSL routines CSINT which evaluate the coefficients of cubic spline and CSDER which calculate the derivative. Further, I calculated the first derivative by using the last suggestion of Jim (by calculating slope at two points, then interpolating). The first two calculations gave the more accurate results than the last one.

Many thanks for all of your valuable suggestions.

Cheers, Nazmul

[Anthony_Richards]

...you are welcome. Glad to have been of some assistance.

[NotThatItMatters]

Not to hijack this post, but has anybody coded the equivalent for equally-spaced data referenced in the 1976 paper from the Journal of Applied Meteorology, "Computing Derivatives from Equally Spaced Data" by Ruthroff and Bodtmann?  The paper proposes taking the Fourier sin transform of the de-trended data and then truncating the result after a few terms.  I have been unsuccessfully attempting to implement that logic in some equally-spaced tabular data I have.

[Anthony_Richards]

By 'detrending' I assume you mean you first fit a some function (or functions) to/through the data then subtract the fitted function from the data? This would then leave you with residuals ('noise') which you then smooth by dropping higher frequency terms? Sounds straightforward to me.

[jimdempseyatthecove]

Or... assuming the frequency is unknown, but known to be fixed, you could do something like a Newton's method to determine the time offset and the frequency. Where the first approximation of the frequency can be determined where the zero point crossings are located and then what the separations are. Keep in mind that you may want to test for harmonics.

Then using this first level approximation, and a derivation of Newton's method, search for a minimal difference between method calculated offset and frequency and the data points (use least squares difference). Note, the derivation of Newton's method must manipulate two input parameters as opposed to one, which is normally done for something like square root.

As to if this is more efficient, you'd have to run the code.

Jim Dempsey

[NotThatItMatters]

Actually, aside from the nuts and bolts, the only problem I am having following the paper is a statement in the paper which is attached.  The constant term in the derivative represents the "de-trending," which is straightforward.  The second equation, (14b), which is the discretized derivative seems incorrect in my view.  Should not the derivative involve the derivative of the filtered function g_F?  This is not discussed, and the results I am getting seem to require it, although I must say at present my derivatives are mis-scaled by the time step.  I probably have missed that factor somewhere.

[Anthony_Richards]

It could be just a typo omitting an apostrophe.

[John_Campbell]

Fitting of a SIN function is a key aspect of tide prediction by doing a least squares fit for phase and amplitude. For tides, the approach is made easier by recognising that the set of possible frequencies is dominated by movement of the earth, moon, sun and other planets. This gives a set of known frequencies that dominate the tidal response and so the least squares analysis only has to estimate the amplitude and phase angle associated with each frequency.

Where the frequency is not known, this is best solved by an iterative solution to a non-linear least squares that separately solves for frequency then phase angle.

Iterative least squares can also be very effective for complex functions, when the function derivative is known. It is just a matter of identifying the most effective way to separate the function into linear parts.

The use of the Lagrange 3-point interpolation formula is best applied in the range (x0-h/2 : x0+h/2), while a 4-point formula is suited to (x0:x1). I mentioned the Lagrange shape function approach developed in Finite Element formulations, as it shows a simple way of developing the interpolation formula as a seperate formula for each sample point. If the point spacing is not uniform, the method can still be applied, with some increase in error estimate due to the estimation of dp/dx (termed the Jacobian matrix).  

Back to the original post; There has been no indication of the accuracy required. Even a simple first difference approach could be used to estimate the derivative in the range (x0:x1), applying this to (x0+x1)/2 and then using a linear interpolation of the derivative from this new difference table.

[Arjen_Markus]

Re equation 14b: it definitely seems to be the omission of a quote/prime.

As for determining the cut-off frequency: in the article the authors use a graphical method to determine where the signal ends (in the frequency domain) and where the noise is dominant. I suppose this could be automated, if the spectrum is not too messy.