在两个已知点之间给定弧长，确定抛物线。

Question

让（0,0）和（Xo，Yo）成为直角坐标系上的两个点。我们想确定通过这两个点并且具有一定弧长S的抛物线曲线Y = AX^2 + BX + C。显然，S > sqrt(Xo^2 + Yo^2)。由于曲线必须从（0,0）经过，因此C = 0。因此，曲线方程简化为：Y = AX^2 + BX。如何确定{A,B}，已知{Xo,Yo,S}？有两个解，我想要A>0的解。

我有一个分析解（复杂），它给出了给定{A,B,Xo,Yo}的S值，但在这里问题是反向的...我可以通过数值求解复杂的方程组来继续进行...但也许有一个数值例程能够完全做到这一点？

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

非常感谢 :-)

- pa23057

---

@ReblochonMasque 不一定。任何经过原点的多项式都有 x - 0 或者简单地说是 x 作为根。因此它不能有与 x 无关的常数项。 - EvilTak

通过“S”的弧长，我理解您的意思是由“0, 0”和“x0, y0”形成的弧的长度为“S”？ - EvilTak

这对我来说更像是一个数学问题而不是编程问题。先试着在纸上解出解析解吧？ - Tony Beta Lambda

你读过https://en.wikipedia.org/wiki/Parabola#Length_of_an_arc_of_a_parabola吗？在那个部分的术语中，您知道*s*和*p*并想要计算*f*，作为第一步。对我来说，那里的方程似乎相当超越，因此我不希望得到精确公式。不同的数值方法可能值得尝试。 - MvG

如果你有(A, B, x0, y0) -> S0，并且很难反转，那么这是一个简单的二维最小化问题，你想要最小化abs(S - S0)，所以你可以尝试一个简单的最小二乘算法...实际上它是一维的，因为你有通过点的限制。 - mikuszefski

一个高级的“曲线拟合”库可以做到这一点：从一个家族（比如二次多项式）中找到通过两个给定点并具有给定长度的曲线。当然，这是数学的，不是所有东西;-) - pa23057

Answer 1

请注意，二次方程 a*x0^2 + b*x0 的弧长（线积分）由从 x = 0 到 x = x0 的 sqrt(1 + (2ax + b)^2) 的积分给出。解决积分后，积分值可表示为 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的函数，我们只需要找到s（a）= S时a的值。这就是我最喜欢的寻根算法牛顿-拉夫逊方法再次发挥作用的地方。

使用牛顿-拉夫逊方法查找根的工作算法如下：

对于要获取其根的函数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)。

唯一的问题是计算一个好的初始猜测。由于s(a)的图形的性质，该函数是牛顿-拉弗森方法的一个很好的候选，对于容差/epsilon为1e-10，初始猜测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))

在CodeSkulptor上运行示例。

请注意，由于对函数输入的弧长进行有理逼近，所得到的系数可能会与预期值略有不同。

- EvilTak