Решение линейного уравнения
https://stackoverflow.com/questions/769
-
08-06-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Мне нужно программно решить систему линейных уравнений на 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.
Решение
Правило КрамераиГауссово исключение— два хороших алгоритма общего назначения (см. также Одновременные линейные уравнения).Если вы ищете код, проверьте ГиНаК, Максима, и СимволическийC++ (конечно, в зависимости от ваших лицензионных требований).
РЕДАКТИРОВАТЬ:Я знаю, что ты работаешь в стране C, но я также должен замолвить словечко за СимПи (система компьютерной алгебры на 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 было бы неплохо разработать свои собственные алгоритмы.
Однако вам, возможно, придется беспокоиться о точности, делении на ноль или очень маленьких числах, а также о том, что делать с бесконечным множеством решений.Я предлагаю использовать стандартный пакет числовой линейной алгебры, такой как ЛАПАК.
Взгляните на Фонд Microsoft Solver.
С его помощью вы можете написать такой код:
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.SimpleSolver
Директивы:
Microsoft.SolverFoundation.Services.Directive
Алгоритм:Первобытный
Арифметика:Гибридный
Цена (точная):По умолчанию
Цена (двойная):Крутой край
Основа:Слабый
Число поворотов:3
===Подробности решения===
Цели:
Решения:
а:0.0785250000000004
б:-0,180612500000001
с:-41,6375875
Вы ищете пакет программного обеспечения, который будет выполнять всю работу или фактически выполнять матричные операции и тому подобное и выполнять каждый шаг?
Первый, мой коллега только что использовал Окамл ГЛПК.Это всего лишь оболочка для ГЛПК, но при этом удаляются многие этапы настройки.Однако, похоже, вам придется придерживаться GLPK на языке C.Что касается последнего, спасибо Delicious за сохранение старой статьи, которую я когда-то изучал LP, PDF.Если вам нужна конкретная помощь в дальнейшей настройке, дайте нам знать, и я уверен, что я или кто-нибудь вернется и поможет, но я думаю, что отсюда все будет довольно просто.Удачи!
Набор шаблонных числовых инструментов от NIST есть инструменты для этого.
Одним из наиболее надежных способов является использование QR-разложение.
Вот пример оболочки, в которой я могу вызвать 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
Что касается эффективности во время выполнения, другие ответили лучше, чем я.Если у вас всегда будет одинаковое количество уравнений в качестве переменных, мне нравится Правило Крамера поскольку это легко реализовать.Просто напишите функцию для вычисления определителя матрицы (или используйте уже написанную, я уверен, вы можете найти ее там) и разделите определители двух матриц.
Лично я неравнодушен к алгоритмам Численные рецепты.(Мне нравится версия C++.)
Эта книга расскажет вам, почему алгоритмы работают, а также покажет вам некоторые довольно хорошо отлаженные реализации этих алгоритмов.
Конечно, вы могли бы просто слепо использовать КЛАПАК (Я использовал его с большим успехом), но сначала я бы вручную напечатал алгоритм исключения Гаусса, чтобы хотя бы иметь представление о том, какая работа была проделана для обеспечения стабильности этих алгоритмов.
Позже, если вы будете заниматься более интересной линейной алгеброй, просматривая исходный код Октава ответит на многие вопросы.
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
