【问题标题】：How to Integrate Arc Lengths using python, numpy, and scipy?如何使用 python、numpy 和 scipy 集成弧长？
【发布时间】：2017-07-07 04:54:28
【问题描述】：

在另一个thread 上，我看到有人设法使用mathematica 积分弧的长度。他们写道：

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

如何使用 numpy 和 scipy 的导入将其转移到 python？特别是，我在他的代码中使用“NIntegrate”函数卡在第 4 行。谢谢您的帮助！

另外，如果我已经有了弧长和垂直轴长度，我将如何反转程序以从已知值中吐出原始参数？谢谢！

【问题讨论】：

标签： python numpy scipy automatic-ref-counting ellipse

【解决方案1】：

据我所知scipy 无法执行符号计算（例如符号微分）。您可能想查看http://www.sympy.org 的符号计算包。因此，在下面的示例中，我分析计算导数（Dx(t) 和 Dy(t) 函数）。

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

编辑：问题的第二部分 - 反转问题

根据您知道积分（弧）的值的事实，您现在可以求解确定弧（半轴、角度等）的一个参数。假设您想解决角度。然后您可以使用scipy 中的非线性求解器之一来还原方程quad(theta) - arcval == 0。你可以这样做：

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

然后：

>>> arc(50)
2.531419526553662
>>> invert(arc(50))
50.000031008458365

【讨论】：

非常感谢您的回复！很抱歉，最后我有一个关于四边形函数的快速问题。如果我已经有 2.5314327 的弧长，我将如何在代码中反转该函数？谢谢！
@FairlyFactual 你说的“倒车”是什么意思？这是一个定积分。我猜你想解方程quad(of something) == 2.53...。问题是：你想解决什么问题？
@FairlyFactual 我将尝试以我理解的方式解决您的“反向问题”添加到我的答案中。

【解决方案2】：

如果您更喜欢纯数值方法，可以使用以下准系统解决方案。这对我很有效，因为我有两个输入 numpy.ndarrays、x 和 y，没有可用的函数形式。

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

注意：我以这种方式分隔返回变量只是为了使其更易读和更清楚地理解正在发生的操作。

顺便提一下，曲线的弧长由下式给出：

【讨论】：