弧線長度使用scipy錯誤的結果

我從一個簡單的例子開始，一個半徑爲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)

有什麼建議嗎？

來源

2017-09-23 Jonathan Nathan

我的記憶可能會誤導我，但是你不需要在你的np.sqrt（）函數中使用+1嗎？那就是：

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

來源

2017-09-23 16:43:39

'1'用於功能定義，其中'y'是'x'的函數，而'y'的導數相對於'x'被使用。這種情況是不同的：參數方程中'x'和'y'都是't'的函數，所以在這裏不應該使用'1'。 –

你是對的。因爲我不得不做任何與弧長有關的事情，所以它太長了。我不知道他的代碼可能會遇到什麼問題。 –

在SciPy的導數函數似乎並不與np.cos（）

當測試derivative(np.cos, np.pi/2)返回0.84147098480789639我的電腦上（實際值應爲-1）很好地工作，但我不知道爲什麼

來源

2017-09-23 17:01:50 leyanpan

爲derivative文檔字符串表示"use a central difference formula with spacing DX "，併爲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)

來源

2017-09-23 17:40:10

如果你正確地執行衍生品，你會得到預期的結果：

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)，則代碼會得出正確的結果

來源

2017-09-23 17:48:19 Nathan

您正在將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。

來源

2017-09-23 23:14:32

弧線長度使用scipy錯誤的結果

回答

相關問題