求解线性方程
Posted
技术标签:
【中文标题】求解线性方程【英文标题】：Solving a linear equation 【发布时间】：2010-09-05 06:32:10 【问题描述】：我需要用 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 的最佳近似值。
【问题讨论】：
其他人已经回答了这个问题，但请查看 Kincaid 和 Cheney 所著的数值分析：科学计算的数学一书。这本书主要是关于求解不同的方程组。 【参考方案1】：Cramer's Rule 和 Gaussian Elimination 是两个很好的通用算法（另见Simultaneous Linear Equations）。如果您正在寻找代码，请查看 GiNaC、Maxima 和 SymbolicC++（当然取决于您的许可要求）。
编辑：我知道你在 C 领域工作，但我也必须为SymPy（Python 中的计算机代数系统）说一句好话。你可以从它的算法中学到很多东西（如果你能读懂一点 python）。此外，它在新的 BSD 许可下，而大多数免费数学包都是 GPL。
【讨论】：
实际上，克莱默法则和高斯消元法在现实世界中都不是很好。两者都没有很好的数值特性，也没有太多用于严肃的应用。尝试 LDU 分解。或者更好，不用担心算法，改用 LAPACK。 对于小于 4 的变量，Cramer 规则最适合编写求解代码 imo【参考方案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)
【讨论】：
【参考方案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);
这是输出： ===求解器基础服务报告=== 日期时间：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
【讨论】：
那么，我们可以从中期待哪些数值稳定性属性？由于这不是开源的，因此应该进行尽职调查，并与 LAPACK 等主流库进行基准测试等。与使用专有解决方案相比，它必须具有一些实质性优势。【参考方案4】：对于 3x3 线性方程组，我想推出自己的算法是可以的。
但是，您可能不得不担心准确性、除以零或非常小的数字以及如何处理无限多的解决方案。我的建议是使用标准的数值线性代数包，例如LAPACK。
【讨论】：
【参考方案5】：您是否正在寻找一个可以完成这项工作或实际执行矩阵运算等并执行每个步骤的软件包？
第一个，我的一个同事刚刚使用了Ocaml GLPK。它只是GLPK 的包装，但它删除了许多设置步骤。不过，看起来您将不得不坚持使用 C 语言中的 GLPK。对于后者，感谢delicious保存了我前段时间学习LP的一篇旧文章，PDF。如果您需要进一步设置的具体帮助，请告诉我们，我敢肯定，我或其他人会回来提供帮助，但是，我认为从这里开始是相当直接的。祝你好运！
【讨论】：
【参考方案6】：来自 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 在库中定义。
【讨论】：
【参考方案7】：在运行时效率方面，其他人的回答比我好。如果你总是有相同数量的方程作为变量，我喜欢Cramer's rule，因为它很容易实现。只需编写一个函数来计算矩阵的行列式（或者使用已经编写的函数，我相信你可以找到一个），然后将两个矩阵的行列式相除。
【讨论】：
【参考方案8】：就个人而言，我偏爱Numerical Recipes 的算法。 （我喜欢 C++ 版本。）
这本书将教你为什么这些算法会起作用，还会向你展示这些算法的一些经过良好调试的实现。
当然，你可以盲目地使用CLAPACK（我使用它取得了很大的成功），但我会先手动输入一个高斯消元算法，至少对这种工作有一个模糊的概念致力于使这些算法稳定。
以后，如果你在做更有趣的线性代数，看看Octave的源代码会回答很多问题。
【讨论】：
【参考方案9】：从您问题的措辞看来，您的方程式似乎比未知数多，并且您希望最大程度地减少不一致。这通常通过线性回归来完成，它可以最小化不一致性的平方和。根据数据的大小，您可以在电子表格或统计包中执行此操作。 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
【讨论】：
【参考方案10】：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
【讨论】：
如果A(n,n) 为零怎么办？【参考方案11】：
对于一般情况，您可以使用 python 和 numpy 进行高斯消除。然后插入值，得到剩余的值。
【讨论】：
以上是关于求解线性方程的主要内容，如果未能解决你的问题，请参考以下文章