Calculate the arclength

I am trying to estimate a good approximation for the arclength of the curve \$ y = \sqrt{1-x^2}\$ for \$x \in [0, 1]\$. I did this using a Bezièr curve, and made it so that the area under the two curves are alike. In order to estimate the error I want to calculate the differences in arclength. The arclength of a quarter circle is \$\frac14\pi\$ so it boils down to estimating

$$\mathrm{Difference}= \frac{L - \frac12\pi}{\frac12\pi}$$

Where \$L\$ is the arclength of my Bezièr curve. The problem, or rather my question is how to find a good method for estimating this arclength.

I tried using recursion and the Pythagorean theorem. I calculate the distance between \$P\$ and \$Q\$. Then I calculate the distances \$PR\$ and \$QR\$, where \$R\$ lies between \$P\$ and \$Q\$. If \$\mathrm{abs}( d(PR) + d(QR) - d(PQ)) < \mathrm{Error}\$ I return the length \$d(PR) + d(QR)\$.

This is painfully slow. Setting digits = 15 takes about 30 seconds on my computer and digits = 16 takes about a minute.

Have I implemented this method incorrectly, or in a highly inefficient way?

```
from decimal import *
import decimal
from math import pi
getcontext().prec = 100

nu = 2-(66-15*pi)**0.5/float(3)

def bezier(t):
f = (1-t)3 + 3(1-t)2t + 3(1-t)t**2*nu
g = 3(1-t)2tnu + 3(1-t)*t2 + t**3
return (f, g)

def arc_length(start_t = 0, stop_t = 1, error = 10**(-6), length = 0):
mean_t = 0.5*(start_t + stop_t)
P = bezier(start_t)
Q = bezier(mean_t)
R = bezier(stop_t)

length_one = decimal.Decimal(pyt(P, Q))
length_two = decimal.Decimal(pyt(Q, R))
if abs(length_one + length_two - length) < error:
return length_one + length_two
else:
return arc_length(start_t, mean_t, error, length_one) + arc_length(mean_t, stop_t, error, length_two)

def pyt(P, Q):
return ((P[0] - Q[0])2 + (P[1] - Q[1])2)**0.5

if __name__ == '__main__':

digits = 5
error = 10**(-digits)
t1 = 0
t2 = 1
length = arc_length(

weird function name pyth

Did you mean pythagoras? Oh, no, you meant distance.

what does bezier do?

What are f,g? It looks like it returns two values based on t, but it is not clear to me what the expressions signify. How did you get those expressions?

It probably derives somehow from f(x) = sqrt(1-x^2) and its derivatives.

Maybe it would help if you could show us what the code looks like for the arclength of the arc formed by f(x) = 1 between -1 and 1 (it should return 2).

arc length algorithm

Also, the arc-length runs the risk of running into a recursion-error. And it's slow due to function calls.

Suggested implementation (typed on phone, untested):

def arclength(f, low, high, error):
    p_low, p_high = f(low), f(high)
    length = 0
    todo = [(low, high, p_low, p_high, distance(p_low, p_high))]

    while todo:
        segment = todo.pop()
        middle = (segment[0] + segment[1])/2
        p_middle = f(middle)
        d_left = distance(p_middle, segment[2])
        d_right = distance(p_middle, segment[3])

        if abs(d_left + d_right - segment[4]) < error:
            length += d_left + d_right
        else:
            todo.append((low, middle, p_low, p_middle, d_left))
            todo.append((middle, high, p_middle, p_high, d_right))

    return length

This has the following invariants on todo

• all tuples are of the form (t0, t1, f(t0), f(t1), distance(f(t0), f(t1))).

• at each iteration, a segment is split in two

• if the change in value from splitting up is small enough, the resulting length is added to length

• else: the separate segments are both added back to todo.

(The f is in this case bezier, this allows you to replace it with other curve-defining functions like (cos(t), sin(t)))

StackExchange Code Review Q#128825, answer score: 2