Мне нужно программно решить систему линейных уравнений в C, Object 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.
Правило Крамера а также Исключение Гаусса являются двумя хорошими универсальными алгоритмами (см. также Синхронные линейные уравнения). Если вы ищете код, посмотрите GiNaC, Maxima и SymbolicС++ (в зависимости от ваших требований к лицензированию, конечно).
EDIT: Я знаю, что вы работаете на земле C, но я также должен сказать хорошее слово SymPy (a система компьютерной алгебры в Python). Вы можете многому научиться из своих алгоритмов (если вы можете прочитать немного python). Кроме того, он под новой лицензией BSD, в то время как большинство бесплатных математических пакетов являются GPL.
Вы можете решить это с помощью программы точно так же, как вы ее решаете вручную (с умножением и вычитанием, а затем возвращаете результаты обратно в уравнения). Это довольно стандартная математика среднего уровня.
-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
Они представляют собой одно и то же уравнение, умноженное на два, поэтому вы не можете получить решение от них - умножение первого на два, а затем вычитание оставляет вас с истинным, но бесполезным утверждением:
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х3, я думаю, было бы неплохо развернуть свои собственные алгоритмы.
Однако вам, возможно, придется беспокоиться о точности, делении на нуль или действительно небольших количествах и о том, что делать с бесконечным числом решений. Мое предложение состоит в том, чтобы пойти со стандартным пакетом числовых линейных алгебр, например LAPACK.
Взгляните на 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 ===
Дата и время: 20.04.2009 23:29:55
Название модели: По умолчанию
Требуемые возможности: LP
Решить время (мс): 1027
Общее время (мс): 1414
Решить статус завершения: Оптимально
Выбрано Solver: Microsoft.SolverFoundation.Solvers.SimplexSolver
Директивы:
Microsoft.SolverFoundation.Services.Directive
Алгоритм: Primal
Арифметика: гибрид
Ценообразование (точное): По умолчанию
Цена (двойной): SteepestEdge
Основа: Сладкая
Сводный граф: 3
=== Сведения о решении ===
Цели:
Решения:
a: 0.0785250000000004
b: -0.180612500000001
c: -41.6375875
Вы ищете программный пакет, который будет выполнять работу или выполнять операции с матрицей и т.д. и делать каждый шаг?
Первый, мой коллега только что использовал Ocaml GLPK. Это просто обертка для GLPK, но он удаляет много шагов по настройке вещей вверх. Похоже, вам придется придерживаться GLPK, в C. Для последнего, благодаря восхитительному для сохранения старой статьи, которую я когда-либо изучал LP, PDF. Если вам нужна конкретная помощь в настройке, сообщите нам, и я уверен, что я или кто-то побродит обратно и помогу, но, я думаю, это довольно прямо отсюда. Удачи!
Набор шаблонов Numeric Toolkit от NIST имеет инструменты для этого.
Одним из наиболее надежных способов является использование QR Decomposition.
Вот пример обертки, чтобы я мог называть "GetInverse (A, InvA)" в моем коде, и он поместит инверсию в InvA.
void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
{
QR<double> qr(A);
invA = qr.solve(I);
}
Array2D определен в библиотеке.
Из формулировки вашего вопроса кажется, что у вас больше уравнений, чем неизвестных, и вы хотите минимизировать несоответствия. Обычно это делается с линейной регрессией, которая минимизирует сумму квадратов несоответствий. В зависимости от размера данных вы можете сделать это в электронной таблице или в статистическом пакете. 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
В отношении эффективности во время выполнения другие ответили лучше, чем я. Если вы всегда будете иметь такое же количество уравнений, как и переменные, мне нравится Cramer правило, которое легко реализовать. Просто напишите функцию для вычисления детерминанта матрицы (или используйте написанную, я уверен, вы можете ее найти) и разделите детерминанты двух матриц.
Лично я не согласен с алгоритмами Numericical Recipes. (Я люблю издание на С++.)
В этой книге вы узнаете, почему работают алгоритмы, и покажите вам довольно хорошо отлаженные реализации этих алгоритмов.
Конечно, вы могли бы просто слепо использовать CLAPACK (я использовал его с большим успехом), но я бы сперва -типа алгоритма гауссова алгоритма, чтобы, по крайней мере, иметь слабое представление о той работе, которая позволила сделать эти алгоритмы стабильными.
Позже, если вы делаете более интересную линейную алгебру, оглядываясь на исходный код Octave, ответит на множество вопросов.
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