The reverse of finding the arc length

Question

As we know, finding the arc length of a function $f(x)$ from $x=a$ to $x=b$ is straightforward and can be implemented numerically.

For a particular function $f(x)$ that I have, suppose that I have numerically computed its arclength to be approximately $L$ using Simpson's rule.

Now, for a set of values $\ell_k$ such that $0 \le \ell_k \le L$, I want to find $x_k\in [a,b]$ for which the arc length of $f(x)$ from $x=a$ to $x = x_k$ is $\ell_k$.

I am not sure what's the best way to do this numerically but I followed the suggestion here. That is, if $L(x)$ is the arclength function, $L(x) = \int_a^x \sqrt{1+f'(x)^2} dx$ which yields a differential equation $\frac{dx}{dL} = \frac{1}{\sqrt{1+f'(x)^2}}$.

Solving this differential equation numerically is then not much of a problem. I used the initial condition that the $x-$value corresponding to arclength 0 is $x=a$.

Now, inevitably, numerical errors arise. That is, after solving the differential equation numerically, the $x-value$ corresponding to length $L$ is not exactly $b$. Oftentimes, it goes slightly beyond $b$. I used Matlab's ode45 to solve the differential equation.

I am wondering if I can impose two conditions on this equation, i.e. $x(0) = a$ and $x(L) = b$. If so, how can I solve the ODE numerically? It seems like bvp methods don't apply.

Any suggestions?

Answer 1

The naive approach (which is possibly justified and adequate in this particular case due to how you obtained $L$ in the first place) is this:

In evaluating $L$, you did in fact practically evaluate $\ell(x)$ (such that $L = \ell(b)$). What you are looking for is $x(\ell)$, its inverse. Thus my suggestion would be to either

• store a list of all intermediate $\ell(x)$ and look up the $x$ closest to a given $\ell$, if memory permits; or:
• in computing $L$, check if the current $\ell(x)>\ell_k$ (where $i$ increases which each match), and store the corresponding $x$ at each $k$. This means you only need to store $k_\max$, but either you have to know the $\ell_k$ before knowing $L$, or you need to accumulate twice (once for $L$, and then later for the $\ell_k$).

This may appear brute force, but it has two justifications:

• In general, solving a differential equation will require the same amount of steps as computing $L$ (and thereby the $\ell_k$, as given in the second suggestion above).
• The $\ell_k$ will be of the same order of accuracy as $L$. When evaluating $L$ and $\ell_k$ in two entirely different ways, you could have the akward result of having some $\ell_k > L$. This can't happen if computed as above.

If you feel you need more accurate results, I suggest the following:

Say you find that $\ell(x_i) < \ell_k$ and $\ell(x_{i+1}) > \ell_k$. My first advice is to assume that $\ell(x_i)$ is exact, because if it isn't, everything afterwards including $L$ will be wrong as well. So doubting $\ell(x_i)$ at this particular point isn't a very good model - if you doubt $\ell(x_i)$, rerun with a smaller $x_{i+1} - x_{i}$. So $\ell(x_i)$ and $\ell(x_{i+1})$ are accurate, but you need values in between.

• You can rerun the particular interval $f(x_i)$ to $f(x_{i+1})$ at a finer scale. Do note that by this, you would be taking $L$ and $x(\ell_k)$ from two different approaches, hence, the values are not really comparable, and the arc length at the finer scale might differ significantly from $\ell(x_i) - \ell(x_{i+1})$. However you would not likely not notice in practice.
• Alternatively you could try to obtain the intermediate values directly from the Simpson approximation from which you took all other values. This would mean that all values are consistent, but the inversion of the Simpson parabola doesn't seem straightforward. So you could do it numerically, as proposed above, but instead of refining $f$, you could refine the approximated parabola.

All in all you might want to consider using the trapezoid scheme instead of Simpson (first degree Newton-Cotes) because it is easier to interpolate (the intermediate $x_k$ would (consistently with all other numbers!) lie between $x_i$ and $x_{i+1}$ as $\ell_k$ does between $\ell(x_i)$ and $\ell(x_{i+1})$.

All in all, pending contrary opinions, I believe that this approach is as computationally expensive as solving a whole new ODE, and numerically stable in the sense that it maintains the accuracy you've achieved with computing $L$, and gives you $x$ values that are consistent with $L$.

Answer 2

Note that everything you compute are only approximations: you mentioned that you compute the arc length $L$ via Simpson's rule; similarly you solve the ODE via a numerical method. So the question is not so much about numerical errors, rather about mutual consistency of your numerical methods.

If your function $f$ is piecewise linear (or you replace your function by a piecewise linear approximation), you can derive explicit formulas for the arc length function on each linear interval and hence also explicit formulas for the $x_k$.

I do not believe this will be possible for higher order approximations for $f$, e.g. $B$-splines or arbitrary functions.