2
$\begingroup$

I'm getting inaccurate results when computing (plotting) the value of a function against the length of a line integral, where the path comes from a set of numerically calculated points. (In contrast, the integral seems fine when I have an analytic formula.)

For example: The matrix named points1 (which comes from stitching together a few numerical integrations, from NDSolve), composed of a few hundred real-valued triples, should define a smooth path (x(s),y(s),V(x(s),y(s))). I compute (interpolate) a parameterized line, and I plot the third component against the line integral of the first two components.

line1 = Table[ListInterpolation[points1[[1]][[All,jj]], {0,1}], {jj,1,3}]
ParametricPlot[{NIntegrate[Sqrt[line1[[1]]'[tt]^2+line1[[2]]'[tt]^2],{tt,0,ss}],
  line1[[3]][ss]}, {ss,0,1}]

In a few spots the vertical axis has small notches, i.e., changes of direction up and then down. The notches are small, but still visible. They are definitely spurious. I want to get five or so significant digits, and this is giving me three.

Plot with spurious notches

I think the problem might have to do with ListInterpolation giving changes in direction even when the set of points that it is interpolating is monotonic. I get the same problem whether I choose the method Spline or Hermite, and whichever method I give NIntegrate.

Suggestions, anyone?

EDIT: Added code for one location where the problem shows up.

points = {
{0.6325952779137012, -0.7760085305401847, 54.07201697427571}, 
{0.62105807453548, -0.790281367826272, 53.57216747776533}, 
{0.5941956290244017, -0.8236073812563007, 52.209319984152756}, 
{0.5229935177928883, -0.9126291157398833, 47.34249397123951}, 
{0.3191937129061434, -1.1738845637283815, 26.841435108009442}, 
{0.30562102581445766, -1.1916609806787646, 25.358076396649462}, 
{0.2921897615948844, -1.2093010504444346, 23.900356246809537}, 
{0.2658926549639088, -1.243979893428035, 21.091685216353746}, 
{0.21635546253250398, -1.309812352799613, 16.05638285614679}, 
{0.13289109669708532, -1.422146100879459, 8.757698109321453}, 
{0.025168620642993245, -1.5688078584392655, 2.4156787966727435}, 
{0.020773587161218077, -1.5747852074054476, 2.2408427807381663}, 
{0.016541263234193266, -1.5805365848560005, 2.0788822741133828`}, 
{0.008539545100479554, -1.5913950062364555, 1.7898452227991974}, 
{-0.005774836825683922, -1.61075381814465, 1.328719820988212}, 
{-0.028775790485831575, -1.641584474796247, 0.7368957514090972}, 
{-0.03120754151431262, -1.644815864289341, 0.6849477369657428}, 
{-0.0335531257347027,-1.6479262887687476, 0.6367430278603834}, 
{-0.03799906267998585, -1.6538030044313163, 0.5504755732460467}, 
{-0.045993704234286795, -1.664298846988861, 0.41201122449151484}, 
{-0.058970977912545186, -1.681084608190775, 0.23202329470223404}
};
 line = Table[ ListInterpolation[points[[All, jj]], {0, 1}],  {jj, 1, 3} ];

 ParametricPlot[
   {
     NIntegrate[Sqrt[line[[1]]'[tt]^2 + line[[2]]'[tt]^2], {tt, 0, ss}],
     line[[3]][ss]
  },
       {ss, 0, 1},
   PlotRange -> {{0, All}, {0, All}},
   AspectRatio -> 1/2
       ]
$\endgroup$
  • 2
    $\begingroup$ Without the data, it is difficult to offer advice. To keep the problem manageable, consider providing the portion of points near one of the discontinuities in your plot, perhaps near 2.2. Also, I notice that the plot in your question is not the one produced by your code. Please provide the actual plot. $\endgroup$ – bbgodfrey Mar 1 '15 at 13:06
  • $\begingroup$ I edited a few lines of code into the question, with data points for one of the spurious notches. $\endgroup$ – Richard Tasgal Mar 1 '15 at 15:58
  • $\begingroup$ see here mathematica.stackexchange.com/a/59342/2079 for a method to directly find the running path length of the line. Accumulate[Norm /@ Differences@points] may be better than your NIntegrate approach. (I cant test from here.. ) $\endgroup$ – george2079 Mar 1 '15 at 16:34
7
$\begingroup$

The interpolation overshoots the next point and reverses direction.

ParametricPlot[{line[[1]][tt], line[[2]][tt]}, {tt, 0.2, 0.3}, 
 Epilog -> {Point[points[[All, 1 ;; 2]]]}]

Mathematica graphics

You can reduce the interpolation order to 1 or use a centripetal parametrization parametrizeCurve from J.M.'s answer.

parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] := 
 FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; 
  MatrixQ[pts, NumericQ]

tvals = parametrizeCurve[points];
line = Interpolation[Transpose[{tvals, #}]] & /@ Transpose@points;

s[t_?NumericQ] := 
  NIntegrate[Norm[{line[[1]]'[tt], line[[2]]'[tt]}], 
   Evaluate @ DeleteDuplicates @
     Flatten[{tt, 0, Select[tvals, 0 < # < t &], t}]];

ParametricPlot[
  Evaluate@{s[ss], line[[3]][ss]}, 
  {ss, 0, 1}, 
  PlotRange -> {{0, All}, {0, All}}, 
  AspectRatio -> 1/2, 
  Epilog -> {Point@Table[{s[ss], line[[3]][ss]}, {ss, 0, 1, 0.1}]}
]

Mathematica graphics

Including the interpolation grid points tvals in NIntegrate speeds up the integration. For a really fast implementation use

s = NDSolveValue[{ss'[t] == Norm[{line[[1]]'[t],line[[2]]'[t]}], ss[0]==0}, ss, {t, 0, 1}]

which constructs an InterpolatingFunction for the arc length.

$\endgroup$
  • $\begingroup$ Thanks! After a correction -- integrating over the roots of the squares of the derivatives of the x- and y-components -- it worked. $\endgroup$ – Richard Tasgal Mar 2 '15 at 7:53
  • $\begingroup$ @RichardTasgal You're welcome, and thanks for fixing it. At some point I had done the problem in terms of a single interpolation that returned the coordinates, and at the last minute I switched back to your method of line. In the cut-n-paste editing, I must have pasted over the primes and not noticed. $\endgroup$ – Michael E2 Mar 2 '15 at 11:13

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.