使用 scipy 计算弧长的错误结果答案

【问题标题】：wrong result for arc length using scipy使用 scipy 计算弧长的错误结果
【发布时间】：2017-09-24 07:23:18
【问题描述】：

我正在尝试推导任意曲线的长度。

我从一个简单的例子开始，一个半径为 R 的圆。我得到了错误的结果！

结果似乎与 R 的真实结果不同，这可能会提示问题。

以下代码：

from scipy.integrate import quad
from scipy.misc import derivative
import numpy as np

r = lambda t: 1
x = lambda t: r(t)*np.cos(t)
Dx = lambda t: derivative(x, t)
y = lambda t: r(t)*np.sin(t)
Dy = lambda t: derivative(y, t)

print(quad(lambda t: np.sqrt(Dx(t)**2 + Dy(t)**2), 0, 2*np.pi))

结果

(5.287118128162912, 5.869880279799524e-14)

对于 R=1，它应该是 2*pi = 6.28...

对于 R=5 它是

(26.435590640814564, 2.9349401398997623e-13)

有什么建议吗？

【问题讨论】：

标签： python math scipy

【解决方案1】：

derivative 的文档字符串是 "use a central difference formula with spacingdx"，dx 的默认值是 1，这太大了，无法准确近似函数的导数。例如，尝试dx=1e-8。

使用您的代码，但将 Dx 和 Dy 更改为

In [21]: Dx = lambda t: derivative(x, t, dx=1e-8)

In [22]: Dy = lambda t: derivative(y, t, dx=1e-8)

这是我得到的：

In [23]: print(quad(lambda t: np.sqrt(Dx(t)**2 + Dy(t)**2), 0, 2*np.pi))
(6.283185278344876, 1.7738885483822232e-08)

【讨论】：

【解决方案2】：

如果你正确地实现导数，你会得到预期的结果：

from scipy.integrate import quad
from scipy.misc import derivative
import numpy as np

r = lambda t: 1
x = lambda t: r(t)*np.cos(t)
Dx = lambda t: -r(t)*np.sin(t)
y = lambda t: r(t)*np.sin(t)
Dy = lambda t: r(t)*np.cos(t)

print(quad(lambda t: np.sqrt(Dx(t)**2 + Dy(t)**2), 0, 2*np.pi))

错误来自函数 derivative，它使用中心有限差分公式，如文档 https://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.derivative.html 所述，默认步长为 1，这是非常多的。此步长的通常值应为1e-5 to 1e-8 的顺序。如果你强制使用derivative(x, t, dx=1e-5)，你的代码会产生正确的结果

【讨论】：

【解决方案3】：

我的记忆可能会误导我，但你不需要在你的 np.sqrt() 函数中 +1 吗？那就是：

print(quad(lambda t: np.sqrt(1 + Dx(t)**2 + Dy(t)**2), 0, 2*np.pi))

【讨论】：

1 用于函数定义，其中y 是x 的函数，并且使用y 相对于x 的导数。这种情况不同：x 和y 都是t 的函数的参数方程，所以这里不应该使用1。
你是对的。自从我不得不做任何与弧长相关的事情以来，已经太久了。我不知道他的代码有什么问题。

【解决方案4】：

scipy 中的导数函数似乎不适用于 np.cos()

在我的电脑上测试derivative(np.cos, np.pi / 2)时返回0.84147098480789639（实际值应该是-1），但我不知道为什么

【讨论】：

【解决方案5】：

您将sympy 函数与numpy 函数混合使用，但它们的用途不同。要充分利用 sympy，请使用其版本的 cos、sin、pi、sqrt 等。如果您这样做，sympy 可以在需要进行数字工作之前进行更多分析。毕竟，sympy 代表“符号 Python”。

这是一种纯粹的sympy 方式来获得你想要的东西：

from sympy import symbols, sqrt, pi, cos, sin, Derivative, integrate

r, t, x, y, Dx, Dy = symbols('r, t, x, y, Dx, Dy')
r = 1
x, y = r * cos(t), r * sin(t)
Dx, Dy = Derivative(x, t).doit(), Derivative(y, t).doit()
print(integrate(sqrt(Dx**2 + Dy**2), (t, 0, 2*pi)))

打印出来的结果是

2*pi

这是完全正确的。如果您改用r = 5，则打印输出为10*pi，再次完全正确。

如果您需要一个数字的近似结果，只需在打印前将表达式包装在 float() 中，您将得到 6.283185307179586 用于 r = 1 和 31.41592653589793 用于 r = 5。

【讨论】：