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

Мне нужно программно решить систему линейных уравнений в 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.

3.08.2008 18:14:24
Другие люди ответили на это, но посмотрите книгу « Численный анализ: математика научных вычислений » Кинкейда и Чейни. Книга в основном о решении различных систем уравнений.
Matthew 26.09.2008 19:46:42
10 ОТВЕТОВ
РЕШЕНИЕ

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

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

19
3.08.2008 18:42:54
на самом деле ни правило Крамера, ни устранение Гаусса не очень хороши в реальном мире. ни один из них не имеет хороших числовых свойств, и ни один не используется для серьезных приложений. попробуй факторизацию ЛДУ. или лучше, не беспокойтесь об алгоритме и используйте вместо этого LAPACK.
Peter 12.02.2010 19:11:53
для переменных с числом меньше 4 правило Крамера лучше всего подходит для написания кода для решения imo
Ge Rong 1.09.2017 02:24:18

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

-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)
15
1.09.2013 03:47:38

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

Решения:
a: 0.0785250000000004
b: -0.180612500000001
c: -41.6375875

7
21.04.2009 03:33:42
Итак, какие свойства числовой устойчивости мы можем ожидать от этого? Поскольку это не открытый исходный код, он должен быть обеспечен должной осмотрительностью, сравнением с основными библиотеками, такими как LAPACK, и т. Д. Должно быть существенное преимущество, чтобы перевесить необходимость использования проприетарного решения.
Evgeni Sergeev 26.09.2016 11:20:43

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

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

7
20.12.2016 05:58:05

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

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

3
3.08.2008 18:27:13

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

3
25.08.2008 18:53:32

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

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

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

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

2
25.08.2008 18:22:48

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

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

2
29.06.2017 06:48:07
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
1
30.12.2014 15:24:38
Так что, если A(n,n)это ноль?
Evgeni Sergeev 26.09.2016 11:22:56