Řešení lineární rovnice

hlasů
35

Musím programově řešit systém lineárních rovnic v C, Objective C, nebo (v případě potřeby), C ++.

Zde je příklad z rovnice:

-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

Z toho, rád bych se dostat co nejtěsnější přiblížení k a, ba tx.

Položena 03/08/2008 v 19:14
zdroj uživatelem
V jiných jazycích...                            


10 odpovědí

hlasů
17

Cramerovo pravidlo a Gaussova eliminace jsou dva dobré, univerzální algoritmy (viz také lineárních rovnic ). Pokud hledáte pro kód, podívejte se GiNaC , Maxima a SymbolicC ++ (v závislosti na licenčních požadavků, samozřejmě).

EDIT: Já vím, že pracujete v jazyce C souši, ale také musím dát v dobrém slova pro SymPy (počítačová algebra systém v Pythonu). Můžete se dozvědět mnoho z jeho algoritmů (pokud si můžete přečíst trochu pythonu). Také je podle nového BSD licencí, zatímco většina volných matematických balíčky jsou GPL.

Odpovězeno 03/08/2008 v 19:37
zdroj uživatelem

hlasů
15

Můžete vyřešit s programem přesně stejným způsobem jej vyřešit ručně (s násobení a odčítání, a pak krmení výsledky zpět do rovnice). To je docela standardní matematika středoškolský úrovni.

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

Takže můžete skončit s:

a =   0.0785250
b =  -0.1806125
c = -41.6375875

Pokud připojíte tyto hodnoty zpět do A, B a C, zjistíte, že jsou správné.

Trik je použít jednoduchý 4x3 matrice, která redukuje zase na 3x2 matice, pak 2x1, který je „= n“, kde n je skutečný počet. Jakmile jste, že jste krmit do další matice až získat další hodnotu, pak tyto dvě hodnoty do následujícího maticového až dokud jste vyřešili všechny proměnné.

Znáte N odlišné rovnice, můžete vždy platit pro N proměnných. Říkám odlišná, protože tyto dvě nejsou:

 7a + 2b =  50
14a + 4b = 100

Jsou stejné rovnice vynásobit dvěma, takže se nemůžete dostat řešení z nich - vynásobením první dva odečtením nechává vás s pravda, ale k ničemu prohlášení:

0 = 0 + 0

Jako příklad lze uvést, tady je nějaký C kód, který vyjde na rovnic, které máte uložené ve vaší otázce. Prvních několik nezbytných typy, proměnné, podpůrnou funkci pro tisk rovnici, a na začátku 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;

Dále, snížení tří rovnic s třemi neznámými na dvě rovnice se dvěma neznámými:

    // 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 ("");

Dále, snížení dvou rovnic se dvěma neznámými jednoho rovnice o jedné neznámé:

    // 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 ("");

Nyní, když máme vzorec typu number1 = unknown * number2, můžeme jednoduše pracovat neznámé hodnotu unknown <- number1 / number2. Poté, co jste přišel na tuto hodnotu out, nahradit ji do jedné z rovnic se dvěma neznámými a přijít na druhou hodnotu. Pak nahradit obě tyto (nyní známé) neznámých do jednoho z původních rovnic a nyní mají hodnoty pro všechny tři neznámých:

    // 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;
}

Výstupem tohoto kodexu odpovídá dřívější výpočty v této odpovědi:

         >: -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)
Odpovězeno 26/02/2009 v 11:59
zdroj uživatelem

hlasů
7

Pro 3x3 soustavy lineárních rovnic, myslím, že to bude v pořádku nasadit své vlastní algoritmy.

Nicméně, budete muset starat o přesnost, dělení nulou, nebo opravdu malých množstvích a co dělat o nekonečně mnoho řešení. Můj návrh je jít se standardní numerické lineární algebry balíček jako LAPACK .

Odpovězeno 08/08/2008 v 18:59
zdroj uživatelem

hlasů
6

Podívejte se na Microsoft Řešitel Foundation .

S ním můžete napsat kód, jako je tento:

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

Zde je výstup:
=== Řešitel Foundation Service Report ===
Datetime: 04/20/2009 23:29:55
Název modelu: Standardní
Capabilities požádal: LP
Řešení Čas (ms): 1027
Celkový čas (ms): 1414
Řešení dokončení Status: Optimální
Řešitel Vybrané: Microsoft.SolverFoundation.Solvers.SimplexSolver
směrnic:
Microsoft.SolverFoundation.Services.Directive
Algorithm: Primal
Arithmetic: Hybrid
Stanovení ceny (přesný): Default
Cena (double): SteepestEdge
Základ: Slack
Pivot Count: 3
== = Solution Podrobnosti ===
Góly:

Rozhodnutí:
a: ,0785250000000004
b: -0,180612500000001
c: -41.6375875

Odpovězeno 21/04/2009 v 04:33
zdroj uživatelem

hlasů
3

Template Numerical Toolkit od NIST má nástroje pro to, že.

Jeden z nejspolehlivějších způsobů je použít rozklad QR .

Zde je příklad z obalu tak, že mohu nazývat „GetInverse (A, Inva)“ v mém kódu, a to bude klást inverzní do Inva.

void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
   {
   QR<double> qr(A);  
   invA = qr.solve(I); 
   }

Array2D je definován v knihovně.

Odpovězeno 25/08/2008 v 19:53
zdroj uživatelem

hlasů
3

Hledáte software, který bude dělat svou práci, nebo dokonce dělat maticové operace a tak a dělat každý krok?

The první, spolupracovník dolu právě použil ocaml GLPK . Je to jen obálka pro GLPK , ale odstraňuje mnoho kroků nastavení věci do pořádku. Vypadá to, že budete muset držet s GLPK, v C, ačkoli. Za druhé, díky lahodné pro ukládání starý článek jsem se učit LP chvíli zpět, ve formátu PDF . Pokud potřebujete nastavení další konkrétní pomoc, dejte nám vědět a jsem si jistý, já nebo někdo bude putovat zpět a pomoci, ale myslím, že je to docela rovně vpřed odtud. Hodně štěstí!

Odpovězeno 03/08/2008 v 19:27
zdroj uživatelem

hlasů
2

Z hlediska účinnosti je run-time, jiní odpověděli lépe než já Pokud jste vždycky bude mít stejný počet rovnic jako proměnné, mám rád Cramerovo pravidlo , jak je to snadné implementovat. Stačí napsat funkci vypočítat determinant matice (nebo použijte ten, který je již písemně, jsem si jistý, můžete najít jeden venku), a rozdělit determinanty dvou matic.

Odpovězeno 19/09/2008 v 04:36
zdroj uživatelem

hlasů
2

Ze znění vaši otázku, vypadá to, že budete mít více rovnic než neznámých a chcete minimalizovat rozpory. To se obvykle provádí s lineární regresí, která minimalizuje součet čtverců nesrovnalosti. V závislosti na velikosti dat, můžete tak učinit v tabulce nebo ve statistickém balíku. R je vysoce kvalitní, bez obalu, který dělá lineární regrese mezi spoustou dalších věcí. Tam je hodně lineární regrese (a hodně Gotcha je), ale jak je to jednoduché udělat pro jednoduché případy. Zde je příklad R s použitím dat. Všimněte si, že „tx“ je zachytit do svého modelu.

> 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  
Odpovězeno 17/09/2008 v 15:22
zdroj uživatelem

hlasů
1
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
Odpovězeno 30/12/2014 v 15:24
zdroj uživatelem

hlasů
1

Osobně jsem nakloněný algoritmů numerické recepty . (Jsem rád C ++ vydání.)

Tato kniha vás naučí, proč algoritmy pracují, a ukázat vám některé docela dobře-ladil implementace těchto algoritmů.

Samozřejmě, že můžete jen slepě používat CLAPACK (Použil jsem to s velkým úspěchem), ale já bych nejprve ručně zadat algoritmus Gaussova eliminace alespoň mít slabou představu o druhu práce, která odešla do tvorby těchto algoritmů stabilní.

Později, pokud děláte zajímavější lineární algebry, rozhlíží zdrojového kódu Octave odpoví na spoustu otázek.

Odpovězeno 25/08/2008 v 19:22
zdroj uživatelem

Cookies help us deliver our services. By using our services, you agree to our use of cookies. Learn more