C# 二次方程组求解

Question

如何解这类方程？

a *X + b * Y + c *Z = q
d *X + e * Y + f *Z = w
X *X + Y * Y + Z *Z = z

我们正在寻找 X、Y、Z。如果不是最后一行中的平方，则可以将其作为典型的线性方程求解，例如使用点数值中的线性方程，或编写高斯消元法。
但我该如何解决这个问题呢？另外，您知道 .NET 中有哪些库可以解决该方程吗？

Answer 1

这可以被视为 2 个平面和一个球体的一组方程。该解决方案找到两个平面(一条线)的交点，然后找到该线与球体的交点。

可能有 0、1 或 2 个独特的解决方案。

代码是 C 语言，但我认为 OP 可以很容易地转换为 c#

// Eq1: a *X + b * Y + c *Z = q
// Eq2: d *X + e * Y + f *Z = w
// Eq3: X *X + Y * Y + Z *Z = z
typedef struct {
  double x,y,z,s;
} plane_t;

typedef struct {
  double x,y,z;
} point_t;

int Interection_PlanePlaneSphere(point_t XYZ[2], const plane_t *abc, const plane_t *def, double radius) {
  // Find intersection of 2 planes
  // V = abc cross def
  point_t V;  // This is really 3D vector, not a point
  V.x = abc->y*def->z - abc->z*def->y;
  V.y = abc->z*def->x - abc->x*def->z;
  V.z = abc->x*def->y - abc->y*def->x;
  // printf("V (%12g, %12g, %12g)\n", V.x, V.y, V.z);

  // Assume both planes go through z plane, e.g. z = 0 and V.z != 0
  // Code could be adapted to not depend on this assumption.
  // abc->x*P.x + abc->y*P.y = abc->s
  // def->x*P.x + def->y*P.y = def->s
  double det = abc->x*def->y - abc->y*def->x;

  // if Planes are parallel ...
  // Code could be adapted to deal with special case where planes are coincident.
  if (det == 0.0) return 0;  // 
  point_t P;
  P.x = ( abc->s*def->y - def->s*abc->y)/det;
  P.y = (-abc->s*def->x + def->s*abc->x)/det;
  P.z = 0.0;
  // L(t) = P + V*t = intersection of the 2 planes
  // printf("p (%12g, %12g, %12g)\n", P.x, P.y, P.z);

  if (radius < 0) return 0; // bad sphere
  // Find where L(t) is on the sphere, or |L(t)| = radius^2
  // (P.x - V.x*t)^2 + (P.y - V.y*t)^2 + (P.z - V.z*t)^2 = radius^2
  // (V.x^2 + V.y^2 + V.z^2)*t^2 - 2*(P.x*V.x + P.y*V.y + P.z*V.z) + (P.x^2 + P.y^2 + P.z^2) = radius^2;
  // Solve quadratic
  double a, b, c;
  a = V.x*V.x + V.y*V.y + V.z*V.z;
  b = -2*(P.x*V.x + P.y*V.y + P.z*V.z);
  c = P.x*P.x + P.y*P.y + P.z*P.z - radius*radius;
  // printf("abc (%12g, %12g, %12g)\n", a, b, c);
  det = b*b - 4*a*c;
  if (det < 0) return 0; // no solutions
  det = sqrt(det);
  double t;
  t = (-b + det)/(2*a);
  XYZ[0].x = P.x + t*V.x;
  XYZ[0].y = P.y + t*V.y;
  XYZ[0].z = P.z + t*V.z;
  if (det == 0.0) return 1;
  t = (-b - det)/(2*a);
  XYZ[1].x = P.x + t*V.x;
  XYZ[1].y = P.y + t*V.y;
  XYZ[1].z = P.z + t*V.z;
  return 2;
}

void Test() {
  plane_t abcq = {2, -1,  1, 5};
  plane_t defw = {1,  1, -1, 1};
  double z = 100;
  point_t XYZ[2];
  int result = Interection_PlanePlaneSphere(XYZ, &abcq, &defw, sqrt(z));
  printf("Result %d\n", result);
  int i=0;
  for (i=0; i<result; i++) {
    printf("XYZ[%d] (%12g, %12g, %12g)\n", i, XYZ[i].x, XYZ[i].y, XYZ[i].z);
  }
  // Result 2
  // XYZ[0] (           2,      5.41014,      6.41014)
  // XYZ[1] (           2,     -8.41014,     -7.41014)
}