Calculate the arclength

Question

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)**2*t + 3*(1-t)*t**2*nu
 g = 3*(1-t)**2*t*nu + 3*(1-t)*t**2 + 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(0, 1, error)
 print length
 diff = (length - decimal.Decimal(pi)/decimal.Decimal(2))/(decimal.Decimal(pi)/decimal.Decimal(2))
 print diff

• ferada: You might want to look at this SO post for performance of Decimal.
• N3buchadnezzar: Wow moving decimal out of the recursion made a huge speedimprovement. Who could have guessed?
• user110732: Why do you import decimal twice? I think you dont need the **from decimal import * ** since you already have import decimal.
• Gareth Rees: I don't understand the question. If you want to compute \$\sqrt{1 - x^2}\,ドル why not use sqrt(1 - x**2)? Are you actually trying to compute something else, for example the arc length of the curve \$y=\sqrt{1 - x^2}\$? Why are you using decimal when you only want 15 digits of precision? Please update the question so it is clear what you are trying to compute and why.
• N3buchadnezzar: Updated =) \$ {} \$

Answer 1

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)))