0

我有一个参数化的二维曲线： (x,y) = f(t)

函数 f 是任意的但可微分，因此我可以使用标准公式计算出沿曲线任意点的微分弧长 ds。通过数值积分微分弧长公式，我还可以得到从起点到曲线上任意点的总弧长 S(t)。我可以控制计算的准确性。

我想从曲线的起点找到总弧长 S = D 的点 (x,y)。如果实现是在 python 中，那就更好了。我将多次这样做，它是计算应用程序的一部分，我需要严格控制准确性和对收敛的信心。

我不知道求根是否是最好的方法，但我的问题相当于 g(t) = S(t) - D 的求根问题，因为 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

好的，@HYRY，这是一个很大程度上基于你的代码片段。您给了我成功所需的提示：使用“quad”而不是“quadrature”。所以我至少会投票给你的答案，但我想补充一下。

首先，您的代码运行速度很快，但比我所追求的精度差了大约 5 个位置。我在您的示例中添加了 quadtol 和 opttol 以尝试说明求积和求根精度的相互作用。我还添加了一个基于默认高容差的循环以暴露速度差异。

sin 示例在准确性方面比圆圈更敏感。我还添加了一条参数化曲线，其弧长由超几何函数给出，并且注释掉了“brentq”选项，因为 fsolve 在此示例中失败，并且在任何情况下，甚至 brentq 在此任务中都相等或更好。

“正交”很慢，但表现出预期的行为：求根速度、准确性和成功随正交容差而变化。

相比之下，“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 回答