0

パラメータ化された 2D 曲線があります: (x,y) = f(t)

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

曲線の始点から円弧の長さの合計が S = D である点 (x,y) を見つけたいと考えています。実装が Python である場合はさらに良いでしょう。私はこれを何度も行います。これは、精度の厳密な制御と収束のある程度の信頼が必要な計算アプリケーションの一部です。

ルート検索が最良のアプローチであるかどうかはわかりませんが、私の質問は、関数 g(t) が S( t) そうではありません。不正確な関数評価は、精度だけでなく、g(t) の単調性も台無しにします。最初から密な数値積分を試みましたが、永遠にかかります。必要な許容範囲に収束することはかなり確実です。根探索アルゴリズムは、進行するにつれて積分精度を遅延制御する必要があり、最初はずさんな評価を要求し、根アルゴリズムが収束するにつれて精度が向上します。

そんな簡単に手に入るものってあるの？それを行う別の賢い方法はありますか？

助けに感謝します

4

2 に答える 2

2

コードを投稿して、何が問題なのか教えていただけますか？

これがtを計算する私のバージョンです。ここでS（t）== D：

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)
于 2012-04-26T07:14:57.757 に答える
0

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

まず、コードは高速に実行されましたが、私が求めていた精度よりも約 5 桁足りませんでした。あなたの例に quadtol と opttol を追加して、求積法と求根精度の相互作用を説明しようとしました。また、速度の違いを明らかにするために、デフォルトの高許容値に基づいてループを追加しました。

sin の例は、精度に関して円よりもはるかに敏感です。また、弧の長さが超幾何関数によって与えられるパレマー化された曲線を追加し、「brentq」オプションをコメント アウトしました。

"quadrature" は遅いですが、予想される動作を示します: 求根速度、精度、および成功率は、直交許容値によって変化します。

対照的に、「quad」は要求された許容範囲を無視し、常により正確な回答を生成するようです。この要求されていない正確さは、私の質問がもはや興味深いかどうか確信が持てないほど高速に動作することを除けば、煩わしいか、説明を招くでしょう。ありがとう！

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)
于 2012-04-27T17:11:47.917 に答える