algorithm - 通过代码求解联立方程

Question

这似乎是一个非常简单和愚蠢的问题，但我发现的关于它的一切都太复杂了，我无法理解。

我有这两个非常基本的联立方程:

X = 2x + 2z
Y = z - x

鉴于我知道 X 和 Y，我将如何找到 x 和 z？手工操作非常容易，但我不知道如何用代码完成。

Answer 1

This seems like an incredibly simple and silly question to ask

完全没有。这是一个很好的问题，不幸的是它的答案很复杂。让我们来解决

a * x + b * y = u
c * x + d * y = v

我在这里坚持 2x2 的情况。更复杂的情况将需要您使用库。

首先要注意的是 Cramer formulas不好用。当您计算行列式时

a * d - b * c

只要你有a * d ~ b * c，那么你就有了catastrophic cancellation .这种情况很典型，你必须防范。

简单性/稳定性之间的最佳权衡是 partial pivoting .假设 |a| > |c|。那么系统就相当于

a * c/a * x + bc/a * y = uc/a
      c * x +    d * y = v 

这是

cx + bc/a * y = uc/a
cx +       dy = v  

现在，将第一个减去第二个 yield

cx +       bc/a * y = uc/a
     (d - bc/a) * y = v - uc/a

现在很容易解决:y = (v - uc/a)/(d - bc/a) 和 x = (uc/a - bc/a * y )/c。计算 d - bc/a 比 ad - bc 更稳定，因为我们除以最大数(它不是很明显，但它成立 - 用非常接近的系数，你会明白它为什么有效)。

现在，如果 |c| > |a|，您只需交换行并以类似方式继续。

在代码中(请检查 Python 语法):

def solve(a, b, c, d, u, v):
    if abs(a) > abs(c):
         f = u * c / a
         g = b * c / a
         y = (v - f) / (d - g)
         return ((f - g * y) / c, y)
    else
         f = v * a / c
         g = d * a / c
         x = (u - f) / (b - g)
         return (x, (f - g * x) / a)

您可以使用完全旋转(要求您交换 x 和 y，以便第一个除法始终是最大的系数)，但这写起来更麻烦，并且对于 2x2 的情况几乎不需要。

对于 n x n 的情况，所有旋转的东西都被封装到 LU decomposition 中，你应该为此使用一个库。