确定两个已知点之间给定弧长的抛物线答案

【问题标题】：Determine parabola with given arc length between two known points确定两个已知点之间给定弧长的抛物线
【发布时间】：2018-02-02 12:05:17
【问题描述】：

令 (0,0) 和 (Xo,Yo) 为笛卡尔平面上的两个点。我们要确定抛物线，Y = AX^2 + BX + C，它从这两点经过并且给定的弧长等于 S。显然，S > sqrt(Xo^2 + Yo^2)。由于曲线必须从 (0,0) 通过，它应该是 C=0。因此，曲线方程简化为：Y = AX^2 + BX。知道 {Xo,Yo,S} 如何确定 {A,B}？有两种解决方案，我想要一个 A>0 的解决方案。

我有一个解析解（复数），它为给定的一组 {A,B,Xo,Yo} 给出 S，尽管这里的问题是相反的......我可以通过数值求解一个复杂的方程组来继续。 .. 但也许有一个数字程序可以做到这一点？

任何有用的 Python 库？其他想法？

非常感谢:-)

【问题讨论】：

@ReblochonMasque 不一定。任何通过原点的多项式都以x - 0 或简单地以x 为根。因此它不能有一个独立于x 的常数项。
弧长为S，我假设你的意思是0, 0和x0, y0形成的弧长是S？
对我来说，这似乎更像是一个数学问题而不是编程问题。先在纸上试出解析解？
你读过en.wikipedia.org/wiki/Parabola#Length_of_an_arc_of_a_parabola吗？在该部分的命名法中，您知道 s 和 p 并且想要计算 f，作为第一步。那里的方程式对我来说似乎是超然的，所以我不希望有一个确切的公式。不同的数值方法可能值得尝试。
如果你有 (A, B, x0, y0) -> S0 并且它不容易反转，这是一个简单的 2D 最小化问题，你想最小化 ´ abs( S - S0 )`，所以你可以尝试一个简单的平方算法....实际上它是一维的，因为你有它必须通过的点的限制

标签： python geometry curve-fitting

【解决方案1】：

请注意，二次a*x0^2 + b*x0 的弧长（线积分）由sqrt(1 + (2ax + b)^2) 从x = 0 到x = x0 的积分给出。在求解积分时，积分的值为0.5 * (I(u) - I(l)) / a，其中u = 2ax0 + b； l = b;和I(t) = 0.5 * (t * sqrt(1 + t^2) + log(t + sqrt(1 + t^2))，sqrt(1 + t^2) 的积分。

自从y0 = a * x0^2 + b * x0，b = y0/x0 - a*x0。将b 的值替换为u 和l、u = y0/x0 + a*x0、l = y0/x0 - a*x0。将u和l代入线积分（弧长）的解中，我们得到弧长作为a的函数：

s(a) = 0.5 * (I(y0/x0 + a*x0) - I(y0/x0 - a*x0)) / a

现在我们有了作为a 的函数的弧长，我们只需要找到a 的值，s(a) = S 的值。这就是我最喜欢的寻根算法 Newton-Raphson method 再次发挥作用的地方。

Newton-Raphson 求根方法的工作算法如下：

对于要获取根的函数f(x)，如果x(i)是对根的i的猜测，

x(i+1) = x(i) - f(x(i)) / f'(x(i))

其中f'(x) 是f(x) 的导数。这个过程一直持续到两次连续猜测之间的差异很小。

在我们的例子中，f(a) = s(a) - S 和 f'(a) = s'(a)。通过链式法则和商法则的简单应用，

s'(a) = 0.5 * (a*x0 * (I'(u) + I'(l)) + I(l) - I(u)) / (a^2)

I'(t) = sqrt(1 + t^2).

剩下的唯一问题是计算一个好的初始猜测。由于the graph of s(a) 的性质，该函数是Newton-Raphson 方法的绝佳候选者，对于1e-10 的容差/epsilon，y0 / x0 的初始猜测在大约 5-6 次迭代中收敛于解。

一旦找到a 的值，b 就是简单的y0/x0 - a*x0。

将其放入代码中：

def find_coeff(x0, y0, s0):
    def dI(t):
        return sqrt(1 + t*t)

    def I(t):
        rt = sqrt(1 + t*t)
        return 0.5 * (t * rt + log(t + rt))

    def s(a):
        u = y0/x0 + a*x0
        l = y0/x0 - a*x0
        return 0.5 * (I(u) - I(l)) / a

    def ds(a):
        u = y0/x0 + a*x0
        l = y0/x0 - a*x0
        return 0.5 * (a*x0 * (dI(u) + dI(l)) + I(l) - I(u)) / (a*a)

    N = 1000
    EPSILON = 1e-10
    guess = y0 / x0

    for i in range(N):
        dguess = (s(guess) - s0) / ds(guess)
        guess -= dguess
        if abs(dguess) <= EPSILON:
            print("Break:", abs((s(guess) - s0)))
            break
        print(i+1, ":", guess)

    a = guess
    b = y0/x0 - a*x0

    print(a, b, s(a))

Run the example on CodeSkulptor.

请注意，由于示例中作为函数输入给出的弧长的有理近似值，因此获得的系数可能与预期值略有不同。

【讨论】：