求解线性方程

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.

Answer 1

Cramer's Rule 和 Gaussian Elimination 是两个很好的通用算法(也见同时线性方程).如果您正在寻找代码,请查看GiNaC,Maxima和SymbolicC++(当然,这取决于您的许可要求).

编辑:我知道你在C领域工作,但我也必须为SymPy(Python中的计算机代数系统)提供一个好词.你可以从它的算法中学到很多东西(如果你能读懂一些python).此外,它是在新的BSD许可下,而大多数免费数学包是GPL.

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

它们是相同的等式乘以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)

Answer 3

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

但是,您可能不得不担心准确性,除以零或非常小的数字以及如何处理无限多的解决方案.我的建议是使用标准的数值线性代数包,如LAPACK.

Answer 4

看一下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 Service Report ===
Datetime:04/20/2009 23:29:55
Model Name:Default
Capabilities required :LP
Solve Time(ms):1027
Total Time(ms):1414
Solve完成状态:
选择最佳求解器:Microsoft.SolverFoundation.Solvers.SimplexSolver
指令:
Microsoft.SolverFoundation.Services.Directive
算法:原始
算术:混合
定价(精确):默认
定价(双精度):SteepestEdge
基础:松弛
枢轴数:3
== =解决方案详细信息===
目标:

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