4
\$\begingroup\$

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
\$\endgroup\$
  • 1
    \$\begingroup\$ You might want to look at this SO post for performance of Decimal. \$\endgroup\$ – ferada May 20 '16 at 8:31
  • \$\begingroup\$ Wow moving decimal out of the recursion made a huge speedimprovement. Who could have guessed? \$\endgroup\$ – N3buchadnezzar May 20 '16 at 8:44
  • \$\begingroup\$ Why do you import decimal twice? I think you dont need the **from decimal import * ** since you already have import decimal. \$\endgroup\$ – user110732 May 21 '16 at 15:39
  • \$\begingroup\$ 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. \$\endgroup\$ – Gareth Rees May 30 '16 at 9:41
  • \$\begingroup\$ Updated =) \$ {} \$ \$\endgroup\$ – N3buchadnezzar May 30 '16 at 9:45
2
\$\begingroup\$

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

\$\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.