# Solving an equation involving ArcLength

I need an efficient way to divide the length of a curve into equal curve length intervals (computational chemistry stuff).

I tried this simple code but it doesn't seem to work and just sits there "running":

NS = number of intervals from 0 to 1;
F = [some function];
L = ArcLength[F, {x, 0, 1}];
For[i = 1, i < NS, i++,
Solve[ArcLength[F, {x, 0, l]/L == i/NS, l]]


Perhaps Mathematica is struggling with the combination of differentiation and integration in an equation to solve.

• please provide reproducible code. What is F, exactly? – Rolf Mertig Jun 15 '18 at 16:20
• The nature of the function doesn't matter for the problem at hand. – gsurfer04 Jun 16 '18 at 12:23

Let's start with a simple curve in the plane.

γ = t \[Function] (2 + Cos[10 Pi t])/3 {Cos[2 Pi t], Sin[2 Pi t]};
ParametricPlot[γ[t], {t, 0, 1}] Compute the velocities of the curve.

Ns = 130;
tlist = Subdivide[0., 1., Ns 1000];
Δt = tlist[] - tlist[];
v = t \[Function] Evaluate[N@Sqrt[\[Gamma]'[t].\[Gamma]'[t]]];
vlist = v /@ tlist;


For each element t of tlist, we compute the distance of γ[t] to γ (using Tai's method ;) ) and store the results in slist.

slist = Join[{0.}, Accumulate[MovingAverage[vlist, 2] Δt]];


Now we create a interpolation function that is an approximation to the inverse of the function

$$s \colon {[0,1]} \to {[0, \operatorname{Length}(\gamma)]}, \qquad s(t) = \int_0^t |\gamma'(r)| \, \operatorname{d}\! r.$$

sinverse = Interpolation[Transpose@{slist, tlist}];


This is how it looks:

Plot[sinverse[t], {t, slist[], slist[[-1]]}, AxesLabel -> {"s", "t"}] We use this function to determine the intervals in the $t$-domain:

intervals = Partition[sinverse /@ Subdivide[0., slist[[-1]], Ns], 2, 1];


And here a test:

intervallengths = ArcLength[γ[t], {t, #[], #[]}] & /@ intervals;
L = ArcLength[γ[t], {t, 0., 1.}];
Max[Abs[intervallengths-L/Ns]]


1.55839*10^-9

If you need it really accurate, you can perform several iterations of Newton's method afterwards:

method = Method -> {"NIntegrate",
PrecisionGoal -> MachinePrecision,
WorkingPrecision -> 20
};
L = ArcLength[γ[t], {t, 0., 1.}, method];
F = x \[Function] ArcLength[γ[t], {t, #[], #[]}, method] & /@ Partition[Join[{0.}, x], 2, 1] - L/Ns;
F' = x \[Function] With[{a = v /@ x},
SparseArray[{Band[{1, 1}] -> a,
Band[{2, 1}] -> -Most[a]}, {Length[x], Length[x]}]
];

x = FixedPoint[
# - LinearSolve[F'[#], F[#], Method -> "Banded"] &,
intervals[[1 ;; -2, 2]]
];
intervals = Partition[Join[{0.}, x, {1.}], 2, 1];

intervallengths = ArcLength[γ[t], {t, #[], #[]}, method] & /@ intervals;

Max[Abs[intervallengths - L/Ns]],


1.12254*10^-15