2013-01-16 98 views
1

如何求解這類方程?C#求解二次方程組

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

我們正在尋找X,Y,Z。如果不是最後一行的正方形,則可以將其作爲典型的線性方程來求解,例如使用Dot Numerics的Linear Equations或寫入Gauss Elimination。
但我該如何解決這個問題?另外,你是否知道.NET中的任何庫可以解決這個問題?

+0

難道這些真的二次方程? 'ax^2 + bx + c = 0' –

+0

這不是上面鏈接的重複!您無法像正常的LinearEquation(使用例如DotNumerics)那樣將其解析爲其他線性方程 – rank1

+0

這是當版主不仔細閱讀問題並用錯誤參數關閉問題時發生的情況... – rank1

回答

1

這可以看作是2個平面和球體的一組方程。該解決方案找到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) 
}