python - 确定两个已知点之间具有给定弧长的抛物线

Question

设 (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库？其他想法？

非常感谢 :-)

Answer 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 方法的极好候选，初始猜测为 y0 / x0对于 1e-10 的容差/epsilon，在大约 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 .

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