# Line integrals over a path defined by a set of points (rather than an analytic formula)

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.

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
]

• 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. Mar 1, 2015 at 13:06
• I edited a few lines of code into the question, with data points for one of the spurious notches. Mar 1, 2015 at 15:58
• 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.. ) Mar 1, 2015 at 16:34

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]]]}]


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}]}
]


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.

• Thanks! After a correction -- integrating over the roots of the squares of the derivatives of the x- and y-components -- it worked. Mar 2, 2015 at 7:53
• @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. Mar 2, 2015 at 11:13