首页 文章

求解线性方程

提问于
浏览
35

我需要以编程方式解决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

由此,我想得到 abtx 的最佳近似值 .

10 回答

  • 3

    Cramer's RuleGaussian Elimination是两个很好的通用算法(另见Simultaneous Linear Equations) . 如果您正在寻找代码,请查看GiNaCMaximaSymbolicC++(当然,取决于您的许可要求) .

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

  • 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
    

    举个例子,这里's some C code that works out the simultaneous equations that you'重新放在你的问题中 . 首先是一些必要的类型,变量,用于打印方程的支持函数,以及 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)
    
  • 15

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

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

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

    Here is the output:
    === Solver Foundation Service Report ===
    日期时间:2009年4月20日23:29:55
    型号名称:默认
    要求的能力:LP
    求解时间(ms):1027
    总时间(毫秒):1414
    解决完成状态:最佳
    选择的解算器:Microsoft.SolverFoundation.Solvers.SimplexSolver
    指令:
    Microsoft.SolverFoundation.Services.Directive
    算法:原始
    算术:混合
    定价(准确):默认
    定价(双倍):SteepestEdge
    基础:松弛
    枢轴数:3
    ===解决方案详情===
    目标:

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

  • 7

    您是在寻找能够完成工作还是实际执行矩阵操作的软件包,并且每个步骤都执行?

    第一个,我的一个同事刚刚使用Ocaml GLPK . 它只是GLPK的包装器,但它消除了很多设置步骤 . 看起来你不得不坚持使用C语言中的GLPK . 对于后者,感谢美味的保存旧文章,我曾经学过LP一段时间,PDF . 如果您需要进一步设置特定帮助,请告诉我们,我可以从这里直接进行操作 . 祝好运!

  • 17

    来自NIST的Template Numerical Toolkit有这样做的工具 .

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

    这是一个包装器的例子,所以我可以在我的代码中调用“GetInverse(A,InvA)”,它会将反转INVA .

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

    Array2D在库中定义 .

  • 6

    根据你的问题的措辞,你似乎有更多的方程而不是未知数,你想要最小化不一致性 . 这通常通过线性回归来完成,其最小化不一致的平方和 . 根据数据的大小,您可以在电子表格或统计包中执行此操作 . 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
    
  • 1

    在运行时效率方面,其他人的回答都比I好 . 如果你总是会有与变量相同数量的方程式,我喜欢Cramer's rule,因为它已经写好了,我相信你能找到一个),并划分两个矩阵的决定因素 .

  • 1

    就个人而言,我偏爱Numerical Recipes的算法 . (我很喜欢C版 . )

    本书将教您算法的工作原理,并向您展示这些算法的一些非常好的调试实现 .

    当然,你可以盲目地使用CLAPACK(我已经非常成功地使用了它),但我首先手动输入高斯消除算法,至少对于制作这些算法所做的工作有一个微弱的想法 . 稳定 .

    之后,如果你正在做更有趣的线性代数,那么查看Octave的源代码将回答很多问题 .

  • 2
    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
    

相关问题