Solving a linear equation
我需要以编程方式求解 C、Objective C 或(如果需要)C 中的线性方程组。
这是一个方程的例子：
1
2
3
| -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 的最佳近似值。
克莱默法则
和
高斯消除
是两种很好的通用算法(另见联立线性方程)。如果您正在寻找代码，请查看 GiNaC、Maxima 和 SymbolicC(当然取决于您的许可要求)。
编辑：我知道你在 C 领域工作，但我还必须为 SymPy(Python 中的计算机代数系统)说一句好话。你可以从它的算法中学到很多东西(如果你能读懂一点 python)。此外，它在新的 BSD 许可下，而大多数免费数学包都是 GPL。
您可以使用程序来解决这个问题，就像手动解决它一样(通过乘法和减法，然后将结果反馈到方程中)。这是相当标准的中学水平数学。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| -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) |
所以你最终得到：
1
2
3
| a = 0.0785250
b = -0.1806125
c = -41.6375875 |
如果将这些值重新代入 A、B 和 C，您会发现它们是正确的。
诀窍是使用一个简单的 4x3 矩阵，该矩阵依次减少为一个 3x2 矩阵，然后是一个 2x1，即 "a = n"，n 是一个实际数字。一旦你有了它，你将它输入下一个矩阵以获得另一个值，然后将这两个值输入下一个矩阵，直到你解决了所有变量。
如果你有 N 个不同的方程，你总是可以求解 N 个变量。我说不同是因为这两个不是：
1
2
| 7a + 2b = 50
14a + 4b = 100 |
它们是相同的等式乘以 2，因此您无法从它们中得到解决方案 - 将第一个乘以 2，然后减去，您会得到正确但无用的陈述：
举例来说，这里有一些 C 代码可以计算出您在问题中的联立方程。首先是一些必要的类型、变量、打印方程的支持函数，以及 main:
的开头
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
| #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)\
",
desc, e-r, e-a, e-b, e-c, post);
}
int main (void) {
double a, b, c; |
接下来，将三个有三个未知数的方程化简为两个有两个未知数的方程：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
| // 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 (""); |
接下来，将两个有两个未知数的方程化简为一个有一个未知数的方程：
1
2
3
4
5
6
7
8
9
10
| // 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 计算出未知值。然后，一旦您计算出该值，将其代入具有两个未知数的方程之一，并计算出第二个值。然后将这两个(现在已知的)未知数代入原始方程之一，您现在拥有所有三个未知数的值：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
| // Finally, substitute values back into equations.
a = equ3[0].r / equ3[0].a;
printf ("From (F ), a = %12.8lf (G)\
", a);
b = (equ2[0].r - equ2[0].a * a) / equ2[0].b;
printf ("From (D,G ), b = %12.8lf (H)\
", 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)\
", c);
return 0;
} |
该代码的输出与此答案中的早期计算相匹配：
1
2
3
4
5
6
7
8
9
10
11
12
| : -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) |
看看 Microsoft Solver Foundation。
用它你可以写出这样的代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
| 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
求解时间(毫秒)：1027
总时间(毫秒)：1414
解决完成状态：最佳
选择的求解器：Microsoft.SolverFoundation.Solvers.SimplexSolver
指令：
Microsoft.SolverFoundation.Services.Directive
算法：原始
算术：混合
定价(精确)：默认
定价(双倍)： SteepestEdge
基础：松弛
枢轴计数：3
===解决方案详情===
目标：
决定：
一个：0.0785250000000004
b: -0.180612500000001
c: -41.6375875
对于 3x3 线性方程组，我想推出自己的算法是可以的。
但是，您可能不得不担心准确性、除以零或非常小的数字以及如何处理无限多的解决方案。我的建议是使用标准的数值线性代数包，例如 LAPACK。
在运行时效率方面，其他人的回答比我好。如果你总是有相同数量的方程作为变量，我喜欢克莱默规则，因为它很容易实现。只需编写一个函数来计算矩阵的行列式(或者使用一个已经编写好的函数，我相信你可以找到一个)，然后将两个矩阵的行列式相除。
来自 NIST 的模板数值工具包具有执行此操作的工具。
其中一种更可靠的方法是使用 QR 分解。
这是一个package器的示例，这样我就可以在我的代码中调用 "GetInverse(A, InvA)" 并将其放入 InvA 中。
1
2
3
4
5
| void GetInverse(const Array2Ddouble& A, Array2Ddouble& invA)
{
QRdouble qr(A);
invA = qr.solve(I);
} |
Array2D 在库中定义。
您是否正在寻找能够完成工作或实际执行矩阵运算等并执行每个步骤的软件包？
第一个，我的一个同事刚刚使用了 Ocaml GLPK。它只是 GLPK 的package器，但它删除了许多设置步骤。不过，看起来您将不得不坚持使用 C 语言中的 GLPK。对于后者，感谢delicious保存了我前段时间学习LP的一篇旧文章，PDF。如果您需要进一步设置的具体帮助，请告诉我们，我敢肯定，我或其他人会回来提供帮助，但是，我认为从这里开始是相当直接的。祝你好运！
从您问题的措辞看来，您的方程式似乎比未知数多，并且您希望最大程度地减少不一致。这通常通过线性回归来完成，它可以最小化不一致性的平方和。根据数据的大小，您可以在电子表格或统计包中执行此操作。 R 是一个高质量、免费的包，它可以进行线性回归等。线性回归有很多(和很多陷阱)，但对于简单的情况来说它很简单。这是一个使用您的数据的 R 示例。请注意，"tx"是您模型的截距。
1
2
3
4
5
6
7
8
9
10
11
12
| 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 |
就个人而言，我偏爱数字食谱的算法。 (我喜欢 C 版。)
这本书将教你为什么这些算法会起作用，还会向你展示这些算法的一些经过良好调试的实现。
当然，你可以盲目地使用 CLAPACK(我使用它取得了巨大的成功)，但我会先手动输入一个高斯消除算法，至少对已经完成的工作有一个模糊的概念使这些算法稳定。
以后，如果你在做更有趣的线性代数，看看 Octave 的源代码会回答很多问题。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
| 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(n1)
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 |
对于一般情况，您可以使用 python 和 numpy 进行高斯消除。然后插入值并获取剩余的值。
|
