scipyを使用した円弧長の結果が間違っています

Question

私は任意の曲線の長さを導出しようとしています。

私は、単純な例から、半径Rの円から始めます。私は間違った結果を得ます！

この結果は、Rによる真の結果とは異なるように見えますが、これは問題のヒントを与えるかもしれません。

次のコードがR = 5それを2 * PI = 6.28 ...

であるべきであるR = 1

(5.287118128162912, 5.869880279799524e-14) 

で

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

結果、は

(26.435590640814564, 2.9349401398997623e-13) 

です。

Answer 1

私の記憶が誤解を招くかもしれませんが、np.sqrt（）関数内で+1が必要なのではないですか？つまり：scipyのダウンロードで

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

Answer 2

導関数はnp.cos（）

それは私のコンピュータ上の0.84147098480789639返しderivative(np.cos, np.pi/2)をテストする（実際の値が-1でなければなりません）でうまく動作するようには思えません、

Answer 3

がderivativeのdocstringは"use a central difference formula with spacing DX "が言うと、dxのデフォルト値は、あなたの関数の導関数の正確な近似を期待するのはあまりにも巨大である、1である理由が、私は考えています。たとえば、dx=1e-8を試してください。あなたのコードで

が、Dxと

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

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

に変更Dyとここに私が得るものです：

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

Answer 4

あなたは正しくデリバティブを実装する場合は、あなたが期待する結果を得る：

を

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

間違いは、の機能から発生します。は、ドキュメントhttps://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.derivative.htmlに記載されているように、中央の有限差分式を使用しています。デフォルトのステップサイズは1です。このステップサイズの通常の値は、1e-5 to 1e-8の順である必要があります。これを強制すると、derivative(x, t, dx=1e-5)でコードが正しい結果を返します

Answer 5

関数を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()を入力してください。r = 1の場合は6.283185307179586、r = 5の場合は31.41592653589793となります。