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]
Improve this question
$\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
    Commented Apr 16, 2022 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$
    – Gopal Verma
    Commented Apr 16, 2022 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
    Commented Apr 16, 2022 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$
    – Gopal Verma
    Commented Apr 16, 2022 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$
    – Daniel Lichtblau
    Commented Apr 16, 2022 at 16:18

1 Answer 1

Reset to default
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.)

Improve this answer
$\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$
    – Daniel Lichtblau
    Commented Apr 17, 2022 at 3:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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