Решение линейного уравнения

Вы ищете программный пакет, который будет выполнять эту работу или фактически выполняет матричные операции и тому подобное и выполняет каждый шаг?

Первый, мой коллега, просто использовал Ocaml GLPK. Это просто оболочка для ГЛПК, но она удаляет многие шаги по настройке. Однако похоже, что вам придется придерживаться GLPK на C. Что касается последнего, спасибо delicious за сохранение старой статьи, которую я некоторое время назад изучал LP, PDF. Если вам нужна конкретная помощь в настройке, дайте нам знать, и я уверен, что я или кто-то вернусь и помогу, но, я думаю, отсюда все довольно просто. Удачи!

Правило Крамера и Исключение по Гауссу это два хороших универсальных алгоритма (см. также Одновременные линейные уравнения). Если вы ищете код, обратите внимание на GiNaC, Максима и СимволическийС ++ (конечно, в зависимости от ваших лицензионных требований).

Обновлено: Я знаю, что вы работаете в стране C, но я также должен замолвить слово SymPy (система компьютерной алгебры на Python). Вы можете многому научиться из его алгоритмов (если вы немного умеете читать на Python). Кроме того, он находится под новой лицензией BSD, в то время как большинство бесплатных математических пакетов являются GPL.

Думаю, для системы линейных уравнений 3x3 было бы нормально развернуть свои собственные алгоритмы.

Однако вам, возможно, придется беспокоиться о точности, делении на ноль или действительно малых числах и о том, что делать с бесконечным количеством решений. Я предлагаю использовать стандартный пакет числовой линейной алгебры, такой как ЛАПАК.

Лично я неравнодушен к алгоритмам Числовые рецепты. (Мне нравится редакция C++.)

Эта книга научит вас, почему алгоритмы работают, а также покажет несколько хорошо отлаженных реализаций этих алгоритмов.

Конечно, вы могли бы просто вслепую использовать КЛАПАН (я использовал его с большим успехом), но я бы сначала напечатал алгоритм исключения Гаусса, чтобы хотя бы иметь слабое представление о том, какая работа была проделана для создания этих алгоритмов. стабильный.

Позже, если вы будете заниматься более интересной линейной алгеброй, просмотр исходного кода Октава ответит на множество вопросов.

У Шаблон Numerical 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 определен в библиотеке.

Судя по формулировке вашего вопроса, кажется, что у вас больше уравнений, чем неизвестных, и вы хотите свести к минимуму несоответствия. Обычно это делается с помощью линейной регрессии, которая сводит к минимуму сумму квадратов несоответствий. В зависимости от размера данных вы можете сделать это в электронной таблице или в статистическом пакете. 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  

Что касается эффективности времени выполнения, другие ответили лучше, чем я. Если у вас всегда будет одинаковое количество уравнений в качестве переменных, мне нравится Правило Крамера, поскольку его легко реализовать. Просто напишите функцию для вычисления определителя матрицы (или используйте уже написанную, я уверен, что вы ее найдете там) и разделите определители двух матриц.

Вы можете решить это с помощью программы точно так же, как вы решаете это вручную (с умножением и вычитанием, а затем возвращением результатов в уравнения). Это довольно стандартная математика средней школы.

-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)

Взгляните на 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
Название Модели: Default
Запрашиваемые возможности: LP
Время решения (мс): 1027
Общее время (мс): 1414
Статус завершения решения: Оптимально Выбранный решатель: Microsoft.SolverFoundation.Solvers.SimplexSolver
Директивы:
Microsoft.SolverFoundation.Services.Directive
Алгоритм: Primal
Арифметика: Hybrid
Цена (точная): Default
Цена (двойная): SteepestEdge
Основа: Slack
Счетчик поворота: 3
=== Подробности решения ===
Голы:

Решения:
а: 0,0785250000000004
б: -0.180612500000001
с: -41.6375875

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