求解一个线性方程

我需要以编程方式解决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 的最佳近似值。

0
额外 编辑
意见: 1
其他人已经回答了这个问题,但请查看Kincaid和Cheney编写的“数值分析:科学计算的数学”一书。这本书很大程度上是关于解决不同的方程组。
额外 作者 Matthew,

13 答案

在Reutenauer的“自由李代数”中,第7.6.2节:

A direct bijection between primitive necklaces of length n over F and the set of irreducible polynomials of degree n in F[x] may be described as follows: let K be the field with qn elements; it is a vector space of dimension n over F, so there exists in K an element θ such that the set {θ, θq, ..., θqn-1} is a linear basis of K over F.

With each word w = a0...an-1 of length n on the alphabet F, associate the element β of K given by β = a0θ + a1θq + ... + an-1 θqn-1. It is easily shown that to conjugate words w, w' correspond conjugate elements β, β' in the field extension K/F, and that w \mapsto β is a bijection. Hence, to a primitive conjugation class corresponds a conjugation class of cardinality n in K; to the latter corresponds a unique irreducible polynomial of degree n in F[x]. This gives the desired bijection.

25
额外

See section 38.10 "Generating irreducible polynomials from Lyndon words" of http://www.jjj.de/fxt/#fxtbook

4
额外

我相信这样的双色注射是在英国展出的

S. Golomb。不可约多项式,同步代码,原始项链和 分圆代数。在Proc。 Conf组合数学。和其Appl。,第358- 370,Chapel Hill,1969。北卡罗莱纳州新闻。

但我没有立即访问这篇文章 - 我很确定它在那里。

3
额外
一点额外的证据(仍然不确定): springerlink.com/index/P6X9P2BV73L2X2GG.pdf “由于基数k上的Lyndon单词与Fk上的不可约多项式之间存在双射,因此[10]参考Golomb的论文。
额外 作者 CodingWithoutComments,
我似乎也无法访问它,但至少有一篇论文( Berstel“rel =”nofollow noreferrer“> jstor.org/stable/2001573”> Berstel 和Reutenauer )表明这是一个开放的问题,事实上,我有和他们一样的动机来提出这个问题,所以我想我应该仔细阅读本文。
额外 作者 Qiaochu Yuan,

Golomb发明的通信依赖于q ^ n阶的原始元素a的选择。然后,对每个Lyndon单词L =(l_0,l_1,...,l_ {n-1})分配具有作为根的单元多项式元素a ^ {m(L)},其中m在i = 0,1,...,n-1上的l_i * q ^ i的整数和。

3
额外

您是否正在寻找能够完成这项工作或实际执行矩阵运算的软件包,并执行每一步?

第一个,我的同事刚刚使用 Ocaml GLPK 。它只是 GLPK 的包装,但它删除了很多设置的步骤。尽管如此,看起来你必须坚持使用GLPK。对于后者,感谢美味的保存我以前学习LP的一篇旧文章, PDF 。如果您需要进一步设置特定帮助,请告诉我们,我确信,我或某人会回来并帮助,但我认为这是相当直接的。祝你好运!

0
额外

就个人而言,我偏爱数字食谱的算法。 (我喜欢C ++版本。)

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

当然,您可以盲目地使用 CLAPACK (我已经使用它取得了巨大的成功),但我会首先输入一个高斯消去算法,至少对于使这些算法稳定的工作类型有一个微弱的想法。

后来,如果你正在做更有趣的线性代数,看看 Octave 的源代码会回答很多问题。

0
额外

你可以用一个程序解决这个问题,就像你手工解决这个问题一样(用乘法和减法,然后把结果反馈到方程中)。这是非常标准的中学数学。

-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的公式,所以你不能从它们得到一个解决方案 - 将第一个乘以二,然后减去,留下真实但无用的声明:

0 = 0 + 0

举例来说,下面是一些C代码,它可以计算出您在问题中出现的联立方程式。首先需要一些必要的类型,变量,用于打印公式的支持函数以及 main 的开始:

#include 

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 ("");

Now that we have a formula of the type number1 = unknown * number2, we can simply work out the unknown value with unknown <- number1 / number2. Then, once you've figured that value out, substitute it into one of the equations with two unknowns and work out the second value. Then substitute both those (now-known) unknowns into one of the original equations and you now have the values for all three unknowns:

    // 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)
0
额外

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

In terms of run-time efficiency, others have answered better than I. If you always will have the same number of equations as variables, I like Cramer's rule as it's easy to implement. Just write a function to calculate determinant of a matrix (or use one that's already written, I'm sure you can find one out there), and divide the determinants of two matrices.

0
额外

Cramer's Rule and Gaussian Elimination are two good, general-purpose algorithms (also see Simultaneous Linear Equations). If you're looking for code, check out GiNaC, Maxima, and SymbolicC++ (depending on your licensing requirements, of course).

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

0
额外
实际上,克雷默的规则和高斯消除在现实世界中都不是很好。既不具有良好的数值性质,也不用于严重的应用。尝试LDU分解。或更好,不要担心算法,并使用LAPACK。
额外 作者 Peter,
对于小于4的变量,Cramer's Rule最适合用于编写解决代码的问题
额外 作者 Ge Rong,

Template Numerical Toolkit from NIST has tools for doing that.

更可靠的方法之一是使用 QR分解

下面是一个包装器的例子,以便我可以在代码中调用“GetInverse(A,InvA)”,并将其反转为InvA。

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

Array2D在库中定义。

0
额外
什么是 qr.solve(I)中的 I
额外 作者 Ponkadoodle,

查看 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===
Datetime: 04/20/2009 23:29:55
Model Name: Default
Capabilities requested: LP
Solve Time (ms): 1027
Total Time (ms): 1414
Solve Completion Status: Optimal
Solver Selected: Microsoft.SolverFoundation.Solvers.SimplexSolver
Directives:
Microsoft.SolverFoundation.Services.Directive
Algorithm: Primal
Arithmetic: Hybrid
Pricing (exact): Default
Pricing (double): SteepestEdge
Basis: Slack
Pivot Count: 3
===Solution Details===
Goals:

Decisions:
a: 0.0785250000000004
b: -0.180612500000001
c: -41.6375875

0
额外
那么,我们可以期待什么数值稳定性?由于这不是开源的,它应该有尽职调查,对LAPACK等主流图书馆的基准测试。要获得专利解决方案,必须有一些实质性的优势。
额外 作者 Evgeni Sergeev,
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
0
额外
那么如果 A(n,n)是零呢?
额外 作者 Evgeni Sergeev,