如何计算给定三个点的抛物线的顶点

Question

我有三个 X/Y 点，形成一条抛物线。 我只需要计算经过这三个点的抛物线的顶点是什么。 最好是一种快速的方法，因为我必须做很多这样的计算！

“询问科学家”网站提供了这个答案：

抛物线的一般形式由以下方程给出：A * x^2 + B * x + C = y 其中 A、B 和 C 是任意实数常数。 您有三对点，它们是 (x,y) 有序对。 将每个点的 x 和 y 值代入抛物线方程。 您将得到三个未知数（三个常数）的三个线性方程。 然后，您可以轻松地求解这个由三个方程组组成的 A、B 和 C 值，并且您将得到与 3 个点相交的抛物线方程。 顶点是一阶导数为 0 的地方，一点代数给出： ( -B/2A , C - B^2/4A ) 顶点。

很高兴看到用 C# 或 C++ 执行此计算的实际代码。 有人吗？

Answer 1

谢谢 David，我将您的伪代码转换为以下 C# 代码：

public static void CalcParabolaVertex(int x1, int y1, int x2, int y2, int x3, int y3, out double xv, out double yv)
{
    double denom = (x1 - x2) * (x1 - x3) * (x2 - x3);
    double A     = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom;
    double B     = (x3*x3 * (y1 - y2) + x2*x2 * (y3 - y1) + x1*x1 * (y2 - y3)) / denom;
    double C     = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + x1 * x2 * (x1 - x2) * y3) / denom;

    xv = -B / (2*A);
    yv = C - B*B / (4*A);
}

这就是我想要的。 抛物线顶点的简单计算。 稍后我将处理整数溢出。

Answer 2

这实际上只是一个简单的线性代数问题，因此您可以用符号进行计算。 当您代入三个点的 x 和 y 值时，您将得到三个未知数的三个线性方程。

A x1^2 + B x1 + C = y1
A x2^2 + B x2 + C = y2
A x3^2 + B x3 + C = y3

解决这个问题的直接方法是反转矩阵

x1^2  x1  1
x2^2  x2  1
x3^2  x3  1

并将其乘以向量，

y1
y2
y3

结果是......好吧，并不完全那么简单;-)我在 Mathematica 中做到了，这里是伪代码中的公式

denom = (x1 - x2)(x1 - x3)(x2 - x3)
A = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom
B = (x3^2 * (y1 - y2) + x2^2 * (y3 - y1) + x1^2 * (y2 - y3)) / denom
C = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + x1 * x2 * (x1 - x2) * y3) / denom

： ，如果您想以数字方式进行矩阵数学，您通常会转向线性代数系统（例如 ATLAS，虽然我不确定它是否有 C#/C++ 绑定）。

无论如何，一旦您获得了通过这些公式计算出的 A、B 和 C 的值，您只需将它们代入问题中给出的表达式 -B / 2A 和 C - B^2/4A 来计算顶点的坐标。¹

请注意，如果原始三个点的坐标使 denom 成为一个非常大或非常小的数字，则直接进行计算可能会出现严重的数值误差。 在这种情况下，最好稍微修改一下，以避免除以分母，无论如何它们都会抵消：

denom = (x1 - x2)(x1 - x3)(x2 - x3)
a = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2))
b = (x3^2 * (y1 - y2) + x2^2 * (y3 - y1) + x1^2 * (y2 - y3))
c = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + x1 * x2 * (x1 - x2) * y3)

然后顶点的坐标是 -b / 2a 和 ( c - b^2 / 4a) / denom。 还有各种其他情况可能会受益于这样的“技巧”，例如如果 A 非常大或非常小，或者如果 C 几乎等于 B^2 / 4A 因此它们的差异非常小，但我认为这些情况差异很大，最好针对具体情况的后续问题进行全面讨论。

将所有这些转换为您选择的语言的代码留给读者作为练习。 （对于任何使用标准中缀运算符语法的语言来说，这应该是相当简单的，例如 AZDean 在 C# 中显示。）

---

¹在答案的最初版本中，我认为这是显而易见的，但似乎有很多人喜欢明确提及它。

Answer 3

通过直接代入，您可以得到以下三个方程：

A*x1^2+B*x1+C=y1
A*x2^2+B*x2+C=y2
A*x3^2+B*x3+C=y3

您可以通过注意到这相当于矩阵乘积来解决这个问题：

[x1^2 x1 1] [A]   [y1]
|x2^2 x2 1|*|B| = |y2|
[x3^2 x3 1] [C]   [y3]

因此，您可以通过对矩阵求逆并将逆矩阵与右侧向量相乘来得到 A、B 和 C。

我发现，当我发布这篇文章时，John Rasch 已链接到更深入地实际求解矩阵方程的教程，因此您可以按照这些说明来获得答案。 反转 3x3 矩阵非常容易，所以这应该不会太难。

Answer 4

下面是一个 Fortran 代码，它实现了 @david-z 和 @AZDean 的解决方案：

subroutine parabola_vertex(x1, y1, x2, y2, x3, y3, xv, yv)
real(dp), intent(in) :: x1, y1, x2, y2, x3, y3
real(dp), intent(out) :: xv, yv
real(dp) :: denom, A, B, C
denom = (x1 - x2) * (x1 - x3) * (x2 - x3)
A     = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom
B     = (x3**2 * (y1 - y2) + x2**2 * (y3 - y1) + x1**2 * (y2 - y3)) / denom
C     = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + &
            x1 * x2 * (x1 - x2) * y3) / denom
xv = -B / (2*A)
yv = C - B**2 / (4*A)
end subroutine

Answer 5

我做了一些类似于 @piSHOCK 的答案，也是基于 @AZDean 的代码。 如果您需要大量运行它（或者像我一样在 Matlab 中使用它），这可能是最快的。

我的假设是x1 == -1, x2 == 0, x3 == 1。

a = y2 - ( y1 + y3) / 2   % opposite signal compared to the original definition of A
b = (y3 - y1) / 4         % half of the originally defined B

xExtr = b / a
yExtr = y2 + b * xExtr    % which is equal to y2 + b*b / a

Answer 6

def vertex(x1,x2,x3,y1,y2,y3):
    '''Given three pairs of (x,y) points return the vertex of the
         parabola passing through the points. Vectorized and common expression reduced.'''
    #Define a sequence of sub expressions to reduce redundant flops
    x0 = 1/x2
    x4 = x1 - x2
    x5 = 1/x4
    x6 = x1**2
    x7 = 1/x6
    x8 = x2**2
    x9 = -x7*x8 + 1
    x10 = x0*x1*x5*x9
    x11 = 1/x1
    x12 = x3**2
    x13 = x11*x12
    x14 = 1/(x0*x13 - x0*x3 - x11*x3 + 1)
    x15 = x14*y3
    x16 = x10*x15
    x17 = x0*x5
    x18 = -x13 + x3
    x19 = y2*(x1*x17 + x14*x18*x6*x9/(x4**2*x8))
    x20 = x2*x5
    x21 = x11*x20
    x22 = x14*(-x12*x7 + x18*x21)
    x23 = y1*(-x10*x22 - x21)
    x24 = x16/2 - x19/2 - x23/2
    x25 = -x17*x9 + x7
    x26 = x0*x1*x14*x18*x5
    x27 = 1/(-x15*x25 + y1*(x20*x7 - x22*x25 + x7) + y2*(-x17 + x25*x26))
    x28 = x24*x27
    return x28,x15 + x22*y1 + x24**2*x27 - x26*y2 + x28*(-x16 + x19 + x23)

Answer 7

这闻起来像家庭作业。 “问科学家”是正确的。 假设您的 3 个点是 (x1, y1)、(x2, y2) 和 (x3, y3)。 然后，您将得到三个线性方程：

| M11 M12 M13 |   | A |   | Z1 |
| M21 M22 M23 | * | B | = | Z2 |
| M31 M32 M33 |   | C |   | Z3 |

其中 M₁₁ = x₁²、M₁₂ = x₁，M₁₃ = 1，Z₁ = y₁，对于其他两行类似地使用 (x2, y2) 和 (x3, y3) 代替 (x1, y1)。

求解这个由 3 个方程组组成的方程组将得到 A、B 和 C 的解。

Answer 8

在 https://ideone.com/y0SxKU 运行，

#include <iostream>
using namespace std;
// calculate the vertex of a parabola given three points
// https://stackoverflow.com/q/717762/16582

// @AZDean implementation with given x values

void CalcParabolaVertex(int x1, int y1, int x2, int y2, int x3, int y3, double& xv,  double& yv)
{
    double denom = (x1 - x2) * (x1 - x3) * (x2 - x3);
    double A     = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom;
    double B     = (x3*x3 * (y1 - y2) + x2*x2 * (y3 - y1) + x1*x1 * (y2 - y3)) / denom;
    double C     = (x2 * x3 * (x2 - x3) * y1 + x3 * x1 * (x3 - x1) * y2 + x1 * x2 * (x1 - x2) * y3) / denom;

    xv = -B / (2*A);
    yv = C - B*B / (4*A);
}

// @piSHOCK immplementation assuming regular x values ( wrong!!! )

void CalcParabolaVertex2( int y1, int y2, int y3, double& xv,  double& yv)
{
double d1 = y1 - y2;
double d2 = y1 - y3;

double a =    -d1 + 0.5 * d2;
double b = 2 * d1 - 0.5 * d2;
double c = -y1;

xv =    -0.5      * b / a;
yv = c - 0.25 * b * b / a;  
}

// corrected immplementation assuming regular x values

void CalcParabolaVertex3( int y1, int y2, int y3, double& xv,  double& yv)
{
double d1 = y1 - y2;
double d2 = y1 - y3;

double a = d1 - 0.5 * d2;
double b = -2 * d1 + 0.5 * d2;
double c = y1;

xv =    -0.5      * b / a;
yv = c - 0.25 * b * b / a;  
}


int main() {
    double xv, yv;
    CalcParabolaVertex( 0, 100, 1, 500, 2, 200, xv, yv );
    cout << xv <<" "<< yv << "\n";
    CalcParabolaVertex2( 100, 500, 200, xv, yv );
    cout << xv <<" "<< yv << "\n";
    CalcParabolaVertex3( 100, 500, 200, xv, yv );
    cout << xv <<" "<< yv << "\n";
    return 0;
}

我添加了一些针对负向峰值的单元测试：在 https://ideone.com/WGK90S 上实时运行