2012-04-26 9 views
0

Pythonでパラメータ化曲線に沿って所定の円弧距離の場所を見つける: (X、Y)= F(T)Iがパラメータ2D曲線を有する

関数fは任意であるが微分可能であり、したがって、私ができます標準式を使用して、任意の点での曲線に沿った微分弧長dsを求める。また、弧の長さの微分方程式を数値的に積分することによって、曲線の始めから弧の長さS(t)の合計を得ることができます。私は計算の精度を制御することができます。

私は曲線の始めから総弧長S = Dの点(x、y)を探したいと思います。実装がPythonであればさらに優れています。私はこれを何度もやっていきますが、これは計算アプリケーションの一部であり、精度の厳密な制御と収束の確信が必要です。

根の所見が最善のアプローチであるかどうかわかりませんが、g(t)= S(t)-Dの根の発見問題に相当します。正確にはS(t)はそうではないからである。不正確な関数評価は、精度だけでなくg(t)の単調性も混乱させます。私は最初から厳密な数値積分を試みましたが、それは永遠にかかります。私は、ルート探索アルゴリズムが積分精度を遅れて制御しなければならないという要求された許容誤差に収束することを確信しています。ルートアルゴリズムが収束するにつれて、誤った評価を要求し精度は向上します。

すぐに利用できるものはありますか?それを行うには別の賢い方法がありますか?

+0

ちょうど状況を理解しようとしている:tは正しいですか?あなたの知られているものは、開始時刻、開始位置、終了時刻、終了カーブの長さ(t0、x0、y0、tF、S(tF)= D)です。あなたは、その変位の最終位置(xF、yF)を見つけたいと思っています。 xに明示的な関数としてカーブを書くことができますか?y = h(x)? – fraxel

+0

こんにちはfraxel:tは、カーブをパラメータ化するダミー変数です。私は時間としてそれを考えても害はないと思う。ちなみに、私はHYRYが私のアプローチを正確に示すコードを投稿することによって私を助けたと思います。 –

+0

しかし、あなたはtを取り除き、xで明示的な関数を得ることができますか?つまり、y = h(x)? (おそらくあなたができる)もしそうなら、私はこれを行うためのクールな方法があるかもしれません。 – fraxel

答えて

2

あなたには、いくつかのコードを投稿し、それと間違って何教えていただけます助けに感謝しますか?ここで

は、S(T)== D Tを計算私のバージョンです:

from scipy.integrate import quad 
from scipy.optimize import fsolve 
from math import cos, sin, sqrt, pi 

def circle_diff(t): 
    dx = -sin(t) 
    dy = cos(t) 
    return sqrt(dx*dx+dy*dy) 

def sin_diff(t): 
    dx = 1 
    dy = cos(t) 
    return sqrt(dx*dx+dy*dy) 

def curve_length(t0, S, length): 
    return quad(S, 0, t0)[0] - length 

def solve_t(curve_diff, length):  
    return fsolve(curve_length, 0.0, (curve_diff, length))[0] 

print solve_t(circle_diff, 2*pi) 
print solve_t(sin_diff, 7.640395578) 
0

OKを、@HYRYは、ここでは主にあなたに基づいてコードスニペットです。私は成功するために必要なヒントを教えてくれました。「直角」の代わりに「クワッド」を使用してください。だから私はあなたの答えに少なくとも投票しますが、私はその話に追加したいと思います。

まず、コードは速く走っていましたが、私が過ぎた精度は5箇所ほどありませんでした。あなたの例にquadtolとopttolを追加して、直角位相と根の正確さの正確さの相互作用を説明しようとしました。私はまた、速度差を明らかにするためにデフォルトの高い公差に基づいてループを追加しました。

sinの例は、精度よりもはるかに円よりも敏感です。私はまた、弧の長さが超幾何関数によって与えられている麻痺カーブを追加しました。この例ではfsolveが失敗し、brentqでさえもこのタスクで同等かそれ以上の理由で「brentq」オプションがコメントアウトされました。

"quadrature"は遅いですが、予想される振る舞いを示しています。根の探索速度、精度、成功率は直交トレランスで変化します。

対照的に、「クワッド」は要求された公差を無視し、常により正確な答えを生成するようです。この無作為の正確さは、私の質問がもう面白いかどうかわからない例では非常に速く働くことを除いて、迷惑なものになるか、または説明を招くでしょう。ありがとう!

from scipy.integrate import quad, quadrature 
from scipy.optimize import fsolve, brentq 
from math import cos, sin, sqrt, pi, pow 

def circle_diff(t): 
    dx = -sin(t) 
    dy = cos(t) 
    return sqrt(dx*dx+dy*dy) 

def sin_diff(t): 
    dx = 1 
    dy = cos(t) 
    return sqrt(dx*dx+dy*dy) 

def hypergeom_diff(t): 
    """ y = t^5 x = t^3 """ 
    dx = 3*t*t 
    dy = 5*pow(t,4) 
    return sqrt(dx*dx+dy*dy) 

def curve_length(t0, S, length,quadtol): 
    integral = quad(S, 0, t0,epsabs=quadtol,epsrel=quadtol) 
    #integral = quadrature(S, 0, t0,tol=quadtol,rtol=quadtol, vec_func = False) 
    return integral[0] - length 

def solve_t(curve_diff, length,opttol=1.e-15,quadtol=1e-10): 
    return fsolve(curve_length, 0.0, (curve_diff, length,quadtol), xtol = opttol)[0] 
    #return brentq(curve_length, 0.0, 3.2*pi,(curve_diff, length, quadtol), rtol = opttol) 

for i in range(1000): 
    y = solve_t(circle_diff, 2*pi) 

print 2*pi 
print solve_t(circle_diff, 2*pi) 
print solve_t(sin_diff, 7.640395578) 
print solve_t(circle_diff, 2*pi,opttol=1e-5,quadtol=1e-3) 
print solve_t(sin_diff, 7.640395578,opttol = 1e-12,quadtol=1e-6) 
print "hypergeom" 
print solve_t(hypergeom_diff, 2.0,opttol = 1e-12,quadtol=1e-12) 
print solve_t(hypergeom_diff, 2.0,opttol = 1e-12,quadtol=1e-6)