Решение линейного уравнения
Мне нужно программно решить систему линейных уравнений в 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,
11 ответов
Правило Крамера и исключение Гаусса - два хороших алгоритма общего назначения (также см. Одновременные линейные уравнения). Если вы ищете код, ознакомьтесь с GiNaC, Maxima и SymbolicC++ (конечно, в зависимости от ваших требований к лицензированию).
РЕДАКТИРОВАТЬ: Я знаю, что вы работаете на земле C, но я также должен сказать хорошее слово для SymPy (система компьютерной алгебры в 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.
Взгляните на http://code.msdn.microsoft.com/solverfoundation.
С его помощью вы можете написать код так:
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
Статус завершения решения: оптимальный
Выбранный решатель: Microsoft.SolverFoundation.Solvers.SimplexSolver
Директивы:
Microsoft.SolverFoundation.Services.Directive
Алгоритм: Первичный
Арифметика: Гибрид
Цены (точные): по умолчанию
Цены (двойные): SteepestEdge
Основа: слабина
Pivot Count: 3
=== Детали решения ===
Цели:
Решения:
а: 0,0785250000000004
b: -0.180612500000001
с: -41,6375875
Что касается эффективности во время выполнения, другие ответили лучше, чем я. Если у вас всегда будет одинаковое количество уравнений в качестве переменных, мне нравится правило Крамера, поскольку его легко реализовать. Просто напишите функцию для вычисления определителя матрицы (или используйте уже написанную, я уверен, что вы можете найти ее там) и разделите определители двух матриц.
Вы ищете программный пакет, который будет выполнять работу или фактически выполнять матричные операции и тому подобное, и выполнять каждый шаг?
Во-первых, мой коллега только что использовал Ocaml GLPK. Это просто оболочка для GLPK, но она устраняет множество шагов по настройке. Похоже, что вам придется придерживаться GLPK, в C, хотя. Что касается последнего, спасибо за то, что он сохранил старую статью, которую я раньше изучал LP, PDF. Если вам нужна конкретная помощь в настройке, дайте нам знать, и я уверен, что я или кто-то еще вернусь и помогу, но, думаю, отсюда довольно просто. Удачи!
Шаблон Numeric Toolkit от NIST имеет инструменты для этого.
Одним из наиболее надежных способов является использование QR-разложения.
Вот пример обертки, так что я могу вызвать "GetInverse(A, InvA)" в моем коде, и она вставит обратное в InvA.
void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
{
QR<double> qr(A);
invA = qr.solve(I);
}
Array2D определен в библиотеке.
Лично я неравнодушен к алгоритмам Численных Рецептов. (Мне нравится издание C++.)
Эта книга научит вас, почему работают алгоритмы, а также покажет вам несколько хорошо отлаженных реализаций этих алгоритмов.
Конечно, вы могли бы просто слепо использовать CLAPACK (я использовал его с большим успехом), но сначала я бы набрал вручную алгоритм исключения Гаусса, чтобы хотя бы иметь слабое представление о той работе, которая ушла в создание этих алгоритмов. стабильный.
Позже, если вы будете заниматься более интересной линейной алгеброй, поиск по исходному коду Octave ответит на множество вопросов.
Из формулировки вашего вопроса кажется, что у вас больше уравнений, чем неизвестных, и вы хотите минимизировать несоответствия. Обычно это делается с помощью линейной регрессии, которая минимизирует сумму квадратов несоответствий. В зависимости от размера данных вы можете сделать это в электронной таблице или в статистическом пакете. 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
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
В общих случаях вы можете использовать python вместе с numpy для исключения Гаусса. А затем вставьте значения и получите оставшиеся значения.