3

On another thread, I saw someone manage to integrate the length of a arc using mathematica.They wrote: 

In[1]:= ArcTan[3.05*Tan[5Pi/18]/2.23]
Out[1]= 1.02051
In[2]:= x=3.05 Cos[t];
In[3]:= y=2.23 Sin[t];
In[4]:= NIntegrate[Sqrt[D[x,t]^2+D[y,t]^2],{t,0,1.02051}]
Out[4]= 2.53143

How exactly could this be transferred to python using the imports of numpy and scipy? In particular, I am stuck on line 4 in his code with the "NIntegrate" function. Thanks for the help!

Also, if I already have the arc length and the vertical axis length, how would I be able to reverse the program to spit out the original paremeters from the known values? Thanks!

mrhallak
  • 1,138
  • 1
  • 12
  • 27
Fairly Factual
  • 161
  • 3
  • 13

2 Answers2

4

To my knowledge scipy cannot perform symbolic computations (such as symbolic differentiation). You may want to have a look at http://www.sympy.org for a symbolic computation package. Therefore, in the example below, I compute derivatives analytically (the Dx(t) and Dy(t) functions).

>>> from scipy.integrate import quad
>>> import numpy as np
>>> Dx = lambda t: -3.05 * np.sin(t)
>>> Dy = lambda t: 2.23 * np.cos(t)
>>> quad(lambda t: np.sqrt(Dx(t)**2 + Dy(t)**2), 0, 1.02051)
(2.531432761012828, 2.810454936566873e-14)

EDIT: Second part of the question - inverting the problem

From the fact that you know the value of the integral (arc) you can now solve for one of the parameters that determine the arc (semi-axes, angle, etc.) Let's assume you want to solve for the angle. Then you can use one of the non-linear solvers in scipy, to revert the equation quad(theta) - arcval == 0. You can do it like this:

>>> from scipy.integrate import quad
>>> from scipy.optimize import broyden1
>>> import numpy as np
>>> a = 3.05
>>> b = 2.23
>>> Dx = lambda t: -a * np.sin(t)
>>> Dy = lambda t: b * np.cos(t)
>>> arc = lambda theta: quad(lambda t: np.sqrt(Dx(t)**2 + Dy(t)**2), 0, np.arctan((a / b) * np.tan(np.deg2rad(theta))))[0]
>>> invert = lambda arcval: float(broyden1(lambda x: arc(x) - arcval, np.rad2deg(arcval / np.sqrt((a**2 + b**2) / 2.0))))

Then:

>>> arc(50)
2.531419526553662
>>> invert(arc(50))
50.000031008458365
AGN Gazer
  • 8,025
  • 2
  • 27
  • 45
  • Thank you so much for the response! I'm sorry, but I just have one quick question about the quad function at the end. If I already have the arc length of 2.5314327, how would I be able to reverse that function in code? Thanks! – Fairly Factual Jul 07 '17 at 15:14
  • @FairlyFactual What do you mean by "reversing"? This is a definite integral. I would guess that you want to solve the equation `quad(of something) == 2.53...`. The question is: what do you want to solve for? – AGN Gazer Jul 07 '17 at 17:52
  • @FairlyFactual I added my attempt to solve your "reverse problem" the way I understood it to my answer. – AGN Gazer Jul 10 '17 at 17:36
4

If you prefer a pure numerical approach, you could use the following barebones solution. This worked well for me given that I had two input numpy.ndarrays, x and y with no functional form available.

import numpy as np

def arclength(x, y, a, b):
    """
    Computes the arclength of the given curve
    defined by (x0, y0), (x1, y1) ... (xn, yn)
    over the provided bounds, `a` and `b`.

    Parameters
    ----------
    x: numpy.ndarray
        The array of x values

    y: numpy.ndarray
        The array of y values corresponding to each value of x

    a: int
        The lower limit to integrate from

    b: int
        The upper limit to integrate to

    Returns
    -------
    numpy.float64
        The arclength of the curve

    """
    bounds = (x >= a) & (y <= b)

    return np.trapz(
        np.sqrt(
            1 + np.gradient(y[bounds], x[bounds])
        ) ** 2),
        x[bounds]
    )

Note: I spaced the return variables out that way just to make it more readable and clear to understand the operations taking place.

As an aside, recall that the arc-length of a curve is given by:

Arc-length Equation

yanniskatsaros
  • 338
  • 4
  • 6