如何使用 python、numpy 和 scipy 积分弧长?
How to Integrate Arc Lengths using python, numpy, and scipy?
在另一个 thread 上,我看到有人设法使用 mathematica.They 积分弧长 写道:
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 行。谢谢您的帮助!
此外,如果我已经有了弧长和垂直轴长度,我如何能够反转程序以从已知值中吐出原始参数?谢谢!
据我所知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
如果您更喜欢纯数值方法,则可以使用以下准系统解决方案。这对我来说效果很好,因为我有两个输入 numpy.ndarray
s,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]
)
注意:我将 return 变量隔开,只是为了使其更易读,更清楚地理解正在发生的操作。
顺便说一句,回想一下曲线的弧长由下式给出:
在另一个 thread 上,我看到有人设法使用 mathematica.They 积分弧长 写道:
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 行。谢谢您的帮助!
此外,如果我已经有了弧长和垂直轴长度,我如何能够反转程序以从已知值中吐出原始参数?谢谢!
据我所知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
如果您更喜欢纯数值方法,则可以使用以下准系统解决方案。这对我来说效果很好,因为我有两个输入 numpy.ndarray
s,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]
)
注意:我将 return 变量隔开,只是为了使其更易读,更清楚地理解正在发生的操作。
顺便说一句,回想一下曲线的弧长由下式给出: