解决一个线性方程

Question

我需要使用C、Objective-C或者（如果需要）C++来编写程序解决一组线性方程组。

以下是一个方程组的例子：

-44.3940 = a * 50.0 + b * 37.0 + tx
-45.3049 = a * 43.0 + b * 39.0 + tx
-44.9594 = a * 52.0 + b * 41.0 + tx

根据此，我想获得a、b和tx的最佳近似值。

- Adam Ernst

---

其他人已经回答了这个问题，但是你可以看一下Kincaid和Cheney的书《数值分析：科学计算的数学》。这本书主要讲解如何解决不同的方程组。 - Matthew

Answer 1

克莱姆法则和高斯消元法是两个通用的良好算法（还可以看一下同时线性方程组）。如果您正在寻找代码，请查看GiNaC，Maxima和SymbolicC++（当然要根据您的许可证要求而定）。

编辑：我知道您正在使用 C 语言，但我还要为SymPy说几句话（它是 Python 中的计算机代数系统）。您可以从其算法中学到很多东西（如果您能读一点 Python 的话）。此外，它采用了新BSD许可证，而大多数免费数学包都采用GPL。

- Brian Jorgensen

---

实际上，Cramer法则和高斯消元在现实世界中都不是非常好的选择。它们的数值性质都不太好，并且在严肃的应用中都不怎么被使用。建议尝试LDU分解，或者更好的方法是不用担心算法，改用LAPACK。 - Peter

对于小于4的变量，我认为使用克拉默法则编写求解代码最好。 - Ge Rong

Answer 2

你可以通过编程来解决这个问题，方法和手动计算一样（使用乘法和减法，然后将结果反馈回方程式中）。这是相当标准的中学数学。

-44.3940 = 50a + 37b + c (A)
-45.3049 = 43a + 39b + c (B)
-44.9594 = 52a + 41b + c (C)

(A-B): 0.9109 =  7a -  2b (D)
(B-C): 0.3455 = -9a -  2b (E)

(D-E): 1.2564 = 16a (F)

(F/16):  a = 0.078525 (G)

Feed G into D:
       0.9109 = 7a - 2b
    => 0.9109 = 0.549675 - 2b (substitute a)
    => 0.361225 = -2b (subtract 0.549675 from both sides)
    => -0.1806125 = b (divide both sides by -2) (H)

Feed H/G into A:
       -44.3940 = 50a + 37b + c
    => -44.3940 = 3.92625 - 6.6826625 + c (substitute a/b)
    => -41.6375875 = c (subtract 3.92625 - 6.6826625 from both sides)

所以你最终得到：

a =   0.0785250
b =  -0.1806125
c = -41.6375875

如果你将这些值插回到A、B和C中，你会发现它们是正确的。

诀窍在于使用一个简单的4x3矩阵，逐步缩减为3x2矩阵，然后为2x1矩阵，其中“a = n”，n是一个实际数。一旦你有了这个，你就把它输入到下一个矩阵中得到另一个值，然后把这两个值输入到上一个矩阵中，直到解决所有变量为止。

只要你有N个不同的方程，你就可以解决N个变量。我说不同的，因为这两个方程不是不同的：

 7a + 2b =  50
14a + 4b = 100

它们是同一个方程乘以二，因此您无法从中得出解-将第一个乘以二然后减去会留下真实但无用的陈述：

0 = 0 + 0

---

以下是一些 C 代码的示例，可以解决你在问题中提到的同时方程。首先是必要的类型、变量、支持输出等式的函数，以及 main 的开头：

#include <stdio.h>

typedef struct { double r, a, b, c; } tEquation;
tEquation equ1[] = {
    { -44.3940,  50, 37, 1 },      // -44.3940 = 50a + 37b + c (A)
    { -45.3049,  43, 39, 1 },      // -45.3049 = 43a + 39b + c (B)
    { -44.9594,  52, 41, 1 },      // -44.9594 = 52a + 41b + c (C)
};
tEquation equ2[2], equ3[1];

static void dumpEqu (char *desc, tEquation *e, char *post) {
    printf ("%10s: %12.8lf = %12.8lfa + %12.8lfb + %12.8lfc (%s)\n",
        desc, e->r, e->a, e->b, e->c, post);
}

int main (void) {
    double a, b, c;

接下来，将三元方程组缩减为二元方程组：

    // First step, populate equ2 based on removing c from equ.

    dumpEqu (">", &(equ1[0]), "A");
    dumpEqu (">", &(equ1[1]), "B");
    dumpEqu (">", &(equ1[2]), "C");
    puts ("");

    // A - B
    equ2[0].r = equ1[0].r * equ1[1].c - equ1[1].r * equ1[0].c;
    equ2[0].a = equ1[0].a * equ1[1].c - equ1[1].a * equ1[0].c;
    equ2[0].b = equ1[0].b * equ1[1].c - equ1[1].b * equ1[0].c;
    equ2[0].c = 0;

    // B - C
    equ2[1].r = equ1[1].r * equ1[2].c - equ1[2].r * equ1[1].c;
    equ2[1].a = equ1[1].a * equ1[2].c - equ1[2].a * equ1[1].c;
    equ2[1].b = equ1[1].b * equ1[2].c - equ1[2].b * equ1[1].c;
    equ2[1].c = 0;

    dumpEqu ("A-B", &(equ2[0]), "D");
    dumpEqu ("B-C", &(equ2[1]), "E");
    puts ("");

接下来，将两个包含两个未知数的方程式化简为一个包含一个未知数的方程式:

    // Next step, populate equ3 based on removing b from equ2.

    // D - E
    equ3[0].r = equ2[0].r * equ2[1].b - equ2[1].r * equ2[0].b;
    equ3[0].a = equ2[0].a * equ2[1].b - equ2[1].a * equ2[0].b;
    equ3[0].b = 0;
    equ3[0].c = 0;

    dumpEqu ("D-E", &(equ3[0]), "F");
    puts ("");

现在我们有一个类似于 number1 = unknown * number2 的公式，我们可以用 unknown <- number1 / number2 简单地计算未知值。然后，一旦你找到了这个值，将其代入其中一个含有两个未知数的方程中并计算出第二个值。然后将这些(现在已知的)未知数代入最初的一个方程中，你现在就有了三个未知数的值：

    // Finally, substitute values back into equations.

    a = equ3[0].r / equ3[0].a;
    printf ("From (F    ), a = %12.8lf (G)\n", a);

    b = (equ2[0].r - equ2[0].a * a) / equ2[0].b;
    printf ("From (D,G  ), b = %12.8lf (H)\n", b);

    c = (equ1[0].r - equ1[0].a * a - equ1[0].b * b) / equ1[0].c;
    printf ("From (A,G,H), c = %12.8lf (I)\n", c);

    return 0;
}

这段代码的输出与此答案中之前的计算结果相匹配：

         >: -44.39400000 =  50.00000000a +  37.00000000b +   1.00000000c (A)
         >: -45.30490000 =  43.00000000a +  39.00000000b +   1.00000000c (B)
         >: -44.95940000 =  52.00000000a +  41.00000000b +   1.00000000c (C)

       A-B:   0.91090000 =   7.00000000a +  -2.00000000b +   0.00000000c (D)
       B-C:  -0.34550000 =  -9.00000000a +  -2.00000000b +   0.00000000c (E)

       D-E:  -2.51280000 = -32.00000000a +   0.00000000b +   0.00000000c (F)

From (F    ), a =   0.07852500 (G)
From (D,G  ), b =  -0.18061250 (H)
From (A,G,H), c = -41.63758750 (I)

- Pax

Answer 3

请查看Microsoft Solver Foundation。

使用它，您可以编写如下代码：

  SolverContext context = SolverContext.GetContext();
  Model model = context.CreateModel();

  Decision a = new Decision(Domain.Real, "a");
  Decision b = new Decision(Domain.Real, "b");
  Decision c = new Decision(Domain.Real, "c");
  model.AddDecisions(a,b,c);
  model.AddConstraint("eqA", -44.3940 == 50*a + 37*b + c);
  model.AddConstraint("eqB", -45.3049 == 43*a + 39*b + c);
  model.AddConstraint("eqC", -44.9594 == 52*a + 41*b + c);
  Solution solution = context.Solve();
  string results = solution.GetReport().ToString();
  Console.WriteLine(results); 

以下是输出结果：
===Solver Foundation 服务报告===
日期时间：04/20/2009 23:29:55
模型名称：默认
所需功能：线性规划
求解时间（毫秒）：1027
总时间（毫秒）：1414
求解完成状态：最优
选择的求解器：Microsoft.SolverFoundation.Solvers.SimplexSolver
指令：
Microsoft.SolverFoundation.Services.Directive
算法：原始对偶法
算术运算：混合运算
定价（精确）：默认
定价（双倍）：陡峭边缘
基础：松弛变量
旋转次数：3
===解决方案详细信息===
目标：

决策：
a: 0.0785250000000004
b: -0.180612500000001
c: -41.6375875

- Bobby Ortiz

---

那么，我们可以从这个数值稳定性方面期待什么性质呢？由于这不是开源的，因此应该经过尽职调查，并与主流库（如LAPACK）进行基准测试等。一定有一些实质性的优势来抵消选择专有解决方案的劣势。 - Evgeni Sergeev

Answer 4

对于一个3x3的线性方程组，我猜你可以自己编写算法。

然而，你可能需要担心精确度、除零或非常小的数字以及无限多个解的处理。我的建议是使用标准的数值线性代数包，如LAPACK。

- Pall Melsted

Answer 5

NIST模板数值工具包提供了相关工具。

其中一种更可靠的方法是使用QR分解。

这里有一个示例包装器，这样我就可以在我的代码中调用“GetInverse（A，InvA）”，并将逆矩阵放入InvA中。

void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
   {
   QR<double> qr(A);  
   invA = qr.solve(I); 
   }

Array2D在库中定义。

- Baltimark

Answer 6

就运行时效率而言，其他人的回答比我好。如果您始终具有与变量相同数量的方程式，则我喜欢Cramer's rule，因为它很容易实现。只需编写一个计算矩阵行列式的函数（或使用已经编写好的函数，我相信您可以在外面找到一个），然后将两个矩阵的行列式相除即可。

- Thelema

Answer 7

你是否在寻找一个软件包，可以处理矩阵运算等操作，并且可以逐步执行每个步骤？

对于前者，我的一位同事刚刚使用了 Ocaml GLPK。它只是 GLPK 的封装器，但它消除了许多设置步骤。然而，似乎你需要使用 C 中的 GLPK。对于后者，感谢 delicious 保存了我以前用来学习线性规划的旧文章 PDF。如果您需要进一步的具体帮助，请告诉我们，我相信我或其他人会回来帮忙，但我认为从这里开始非常简单。祝你好运！

- nlucaroni

Answer 8

从您问题的措辞来看，似乎您的方程式比未知数多，您想要最小化不一致性。这通常可以使用线性回归来完成，它会最小化不一致性的平方和。根据数据的大小，您可以在电子表格或统计软件中完成此操作。R是一个高质量且免费的软件包，其中包括线性回归等许多功能。关于线性回归有很多需要注意的地方，但对于简单的情况而言，它是直截了当的。以下是使用您的数据的R示例。请注意，“tx”是您模型的截距。

> y <- c(-44.394, -45.3049, -44.9594)
> a <- c(50.0, 43.0, 52.0)
> b <- c(37.0, 39.0, 41.0)
> regression = lm(y ~ a + b)
> regression

Call:
lm(formula = y ~ a + b)

Coefficients:
(Intercept)            a            b  
  -41.63759      0.07852     -0.18061  

- David Nehme

Answer 9

就我个人而言，我偏爱Numerical Recipes的算法。（我喜欢C++版本。）

这本书将教你为什么算法有效，并向你展示那些算法经过了良好调试的实现。

当然，你也可以盲目地使用CLAPACK（我曾经非常成功地使用过），但我建议你首先手动输入高斯消元算法，至少对这些算法的稳定性有一些模糊的概念。

如果你正在进行更有趣的线性代数工作，查看Octave源代码将会解决很多问题。

- Frank Krueger

Answer 10

function x = LinSolve(A,y)
%
% Recursive Solution of Linear System Ax=y
% matlab equivalent: x = A\y 
% x = n x 1
% A = n x n
% y = n x 1
% Uses stack space extensively. Not efficient.
% C allows recursion, so convert it into C. 
% ----------------------------------------------
n=length(y);
x=zeros(n,1);
if(n>1)
    x(1:n-1,1) = LinSolve( A(1:n-1,1:n-1) - (A(1:n-1,n)*A(n,1:n-1))./A(n,n) , ...
                           y(1:n-1,1) - A(1:n-1,n).*(y(n,1)/A(n,n))); 
    x(n,1) = (y(n,1) - A(n,1:n-1)*x(1:n-1,1))./A(n,n); 
else
    x = y(1,1) / A(1,1);
end

- Bulent S.

---

那么，如果 A(n,n) 是零呢？ - Evgeni Sergeev