解线性方程

Cramer's Rule 和 Gaussian Elimination 是两种好的通用算法(另请参见 Simultaneous Linear Equations )如果您正在寻找代码,请检查 GiNaC , Maxima 和 SymbolicC++ (当然,这取决于您的许可要求)。

编辑:我知道你在C区工作,但我也必须为 SymPy (Python中的计算机代数系统)。你可以从它的算法中学到很多东西(如果你能读一点python的话)。此外,它还拥有新的BSD许可证,而大多数免费的数学软件包都是GPL。

你可以用一个程序来解决这个问题,方法和你用手来解决的方法完全一样(用乘法和减法,然后把结果反馈到等式中)。这是相当标准的中学数学。

-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

它们是 相同的 方程式乘以2,这样你就无法从中得到答案——先乘以2,然后减去,这样你就得到了真实但无用的陈述:

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)

对于一个3x3的线性方程组,我想可以推出自己的算法。

但是,您可能需要担心精度、被零或非常小的数字除以及如何处理无穷多的解。我的建议是使用标准的数字线性代数包,例如 LAPACK .

看看 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); 

以下是输出:
===求解器基础服务报告===
日期:2009年4月20日23:29:55
型号名称:默认
要求的能力:LP
解算时间(ms):1027
总时间(ms):1414
求解完成状态:最优
已选择解算器:microsoft.solverfoundation.solvers.simplexsolver
指令:
Microsoft.SolverFoundation.Services.Directive指令
算法:primal
算术:混合
定价(精确):默认
定价(双倍):最高价
依据:松弛
枢轴计数:3
==解决方案详细信息===
目标:

决定:
答:0.0785250000000004
B:-0.180612500000001
C:-41.6375 875

你是在寻找一个能完成工作的软件包,还是在做矩阵运算之类的工作,并完成每一步?

第一个,我的一个同事刚刚用过 Ocaml GLPK . 它只是一个包装 GLPK 但是它消除了很多设置的步骤。不过,在C语言中,看起来你必须坚持使用glpk。对于后者,感谢美味的保存了一篇我以前学过的老文章。 PDF . 如果您需要进一步的帮助,请告诉我们,我确信,我或其他人会回来帮忙,但是,我认为这是相当直接的。祝你好运!

Template Numerical Toolkit 来自NIST的工具可以做到这一点。

更可靠的方法之一是使用 QR Decomposition .

下面是一个包装器的例子,这样我可以在代码中称为“getInverse(a,inva)”,它将把Inverse放入inva中。

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

array2d在库中定义。

从你的问题的措辞来看,似乎你有比未知更多的方程式,你想要最小化不一致。这通常是通过线性回归来完成的,线性回归可以最小化不一致的平方和。根据数据的大小,您可以在电子表格或统计数据包中执行此操作。R是一个高质量的,免费的包,可以做线性回归,还有很多其他的事情。线性回归有很多种方法(还有很多gotcha方法),但是对于简单的情况,这是很简单的。下面是一个使用您的数据的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  

在运行时效率方面,其他人的回答比我好。如果你总是有相同数量的方程作为变量,我喜欢 Cramer's rule 因为它很容易实现。只需编写一个函数来计算一个矩阵的行列式(或者使用一个已经写过的函数,我相信你可以在那里找到一个),然后将两个矩阵的行列式分开。

就我个人而言,我偏爱 Numerical Recipes . (我喜欢C++版。)

这本书将教你为什么算法工作,并向你展示一些非常好的调试这些算法的实现。

当然,你可以盲目地使用 CLAPACK (我已经成功地使用了它),但是我会首先输入一个高斯消元算法,至少对使这些算法稳定所做的工作有一个模糊的概念。

稍后,如果您正在做更有趣的线性代数,请查看 Octave 会回答很多问题。

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