1
$\begingroup$

I am trying to calculate the derivative of interpolated data but it behaves differently from the analytic solution. Detailed code that I have used follows:

 \[Gamma]w = 70.0*10^-3;(*SR in N/m*)\[Rho] = 1000; c = 3*10^8; g = 9.8;
  we1 = 7*10^-6;(*beam waist*)n = 1.33; P0 = 4.0; rng = 50*10^-6;
 lc = Sqrt[\[Gamma]w/(\[Rho]*g)]; P0 = 4; we1 = 
  7*10^-6; P1 = (P0/(\[Gamma]w*c*Pi))*((n - 1)/(n + 1));
  f[r_] := 7*P1*BesselK[0, r/lc];
  lst2 = Table[{r, f[r]}, {r, we1, rng, rng/1000}]; h2 = 
   Interpolation[lst2, InterpolationOrder -> 1];
   pA = Plot[{h2[r]}, {r, we1, rng}, PlotStyle -> {Blue}, 
   PlotRange -> All, Frame -> True];
  pIn = Plot[{f[r]}, {r, we1, rng}, PlotStyle -> {Red}, 
   PlotRange -> All, Frame -> True];
   Show[{pA, pIn}](*height*);
    A[r_] := h2''[r] + 1/r*h2'[r];(*Interpolation*)
     A2 = D[f[r], {r, 2}] + 1/r*D[f[r], r];(*Analytic*)A3 = D[A2, {r, 
      1}];
     plots = Plot[{A'[r], A3}, {r, we1, rng/1}, PlotRange -> All, 
    Frame -> True, PlotRange -> All, Frame -> True, 
    PlotStyle -> {{Red, Thick}, {Green, Thick}}, 
    GridLines -> {{we1, 2 we1}, {0}}, 
   GridLinesStyle -> Directive[Red, Dashed], 
   PlotLegends -> {"Interpolation (A'[r])", "Analytic (A3)"}]
    [![plots][1]][1]
$\endgroup$
9
  • $\begingroup$ You are using a fairly high interpolation order, which makes the derivatives numerically unstable. I would suggest using a smaller InterpolationOrder. But generally, what you are trying to do is numerically unstable. Better to compute the derivatives straight from the data instead of going through an interpolation. $\endgroup$
    – Roman
    Apr 16 at 10:53
  • $\begingroup$ @Roman, Thanks. I also tried with a small InterpolationOrder. Could you suggest how to compute derivative straight from the data? $\endgroup$ Apr 16 at 11:19
  • $\begingroup$ The easiest way to compute numerical derivatives directly from the data is probably to set InterpolationOrder -> 1 and let Mathematica take care of the details. In this way it will interpolate linearly between data points, thus effectively using the finite-differences method of numerical differentiation. In general, to compute $n^{\text{th}}$ derivatives you should use an interpolation of oder $n$ (but not higher). $\endgroup$
    – Roman
    Apr 16 at 14:44
  • $\begingroup$ Thanks @ Roman, I tried with this code ( lst2 = Table[{r, f[r]}, {r, we1, rng, rng/1000}]; h2 = Interpolation[lst2, InterpolationOrder -> 1] ) but I found huge difference in interpolated data plot and analytic one. $\endgroup$ Apr 16 at 15:44
  • 2
    $\begingroup$ Also there is the recently added Wolfram Function Repository function ListD. It might give reasonable results. And the web page has examples that might give ideas for related methods involving local smoothing. $\endgroup$ Apr 16 at 16:18

1 Answer 1

3
$\begingroup$

InterpolationOrder -> 1 means the first derivative will be discontinuous (piecewise constant) and the second and third derivatives will be zero. You need an interpolation order higher than 3 to get a continuous result.

h2 = Interpolation[lst2, InterpolationOrder -> 7]; 
A[r_] := h2''[r] + 1/r*h2'[r];(*Interpolation*)
A2 = 
 D[f[r], {r, 2}] + 1/r*D[f[r], r];(*Analytic*)
A3 = D[A2, {r, 1}];
plots = Plot[{A'[r], A3}, {r, we1, rng/1}(*,PlotRange->All*), 
  Frame -> True(*,PlotRange->All*), Frame -> True, 
  PlotStyle -> {{Red, Thick}, {Green, Thick}}, 
  GridLines -> {{we1, 2 we1}, {0}}, 
  GridLinesStyle -> Directive[Red, Dashed], 
  PlotLegends -> {"Interpolation (A'[r])", "Analytic (A3)"}]

enter image description here

Still a bit of noise...

P.S.

Using chebInterpolation and chebSeries from FunctionInterpolation over an open interval, we also get a pretty good approximation, which is much less noisy:

h2 = chebInterpolation[{{N@{we1, rng}, 
    chebSeries[f, N@{we1, rng}, 64][[;; -12]]}}]

This is nearly equivalent:

lst2 = Table[{r, f[r]}, {r, 
   Rescale[
    Sin[Pi/2. Range[-64, 64, 2]/64], {-1, 1}, {we1, rng}]}]; 
h2 = Interpolation[lst2, InterpolationOrder -> All]

Note: The function A3 seems to lose about 7 digits of precision due to subtractive cancellation. (Based on evaluating a high-precision A3 on arbitrary precision input and observing the precision loss.)

$\endgroup$
1
  • $\begingroup$ Yeah, that precision loss will be difficult to avoid. Offhand I do not know if smoothing or spectral methods can do better or, if so, by how much. $\endgroup$ Apr 17 at 3:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.