2
$\begingroup$

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.

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

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

enter image description here

Compute the velocities of the curve.

Ns = 130;
tlist = Subdivide[0., 1., Ns 1000];
Δt = tlist[[2]] - tlist[[1]];
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 γ[0] (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[[1]], slist[[-1]]}, AxesLabel -> {"s", "t"}]

enter image description here

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, #[[1]], #[[2]]}] & /@ 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, #[[1]], #[[2]]}, 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, #[[1]], #[[2]]}, method] & /@ intervals;

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

1.12254*10^-15

$\endgroup$

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.