Joachim Mohr   Mathematik Musik

15. Lektion: Der Gauß'sche Algorithmus zur Lösung eines linearen Gleichungssystems

Dazu: Mehrdimensionale Arrays.

Hinweis: Lektion 15 wird für Lektion 16 u.s.w. nicht vorausgesetzt.

Da wir das Verfahren auf dynamische Arrays ausdehnen wollen, fangen wir stets mit dem Index 0 an. Statt die Unbekannten mit x1, x2, x3, ..., xn zu bezeichnen, verwenden wir die Bezeichnung x0, x1, x2, ... , xn-1.

Beispiel 15.1 Es soll ein Gleichungssystem mit 3 Gleichungen und 3 Unbekannten gelöst werden.
koeff  *x0 + koeff  *x1 + koeff * x2 = koeff
     00           01           02           03

koeff  *x0 + koeff  *x1 + koeff * x2 = koeff
     10           11           12           13

koeff  *x0 + koeff  *x1 + koeff * x2 = koeff
     20           21           22           23

zum Beispiel:

   5x - 7x +  x3 = 9
     1    2

  -3x + 2x + 3x  = 4
     1    2    3

   2x + 3x + 4x  = 0
     1    2    3


Dafür schreibt man auch kurz nur die "Koeffizienten":

  5   -7  1  9

 -3    2  3  4

  2    3  4  0


Beim Gauß-Verfahren wird nun folgendermaßen gerechnet:


1. Schritt:

   1. Reihe bleibt

   2. Reihe = zweite Reihe + 3/5*erste  Reihe

   3. Reihe = dritte Reihe - 2/5*erste  Reihe


   Damit sie besser lesbar sind, werden sie im folgenden noch mit einer

   passenden Zahl durchmultipliziert. Ergebnis:


   5      -7       1        9

   0       11     -18      -47

   0       29      18      -18


2. Schritt:

   1. Reihe = erste Reihe  + 11/7*zweite Reihe

   2. Reihe bleibt

   3. Reihe = dritte Reihe -11/29*dritte Reihe


   Ergebnis (nach durchmultiplizieren):

  11      0      -23      -46

  0       11     -18      -47

  0       0       1       233/144


3. Schritt:

   1. Reihe = erste Reihe  -23/1*dritte Reihe

   2. Reihe = zweite Reihe +18/1*dritte Reihe

   3. Reihe = bleibt

  Ergebnis (nach durchmultiplizieren):


  1       0       0       -115/144

  0       1       0       -13/8

  0       0       1       233/144


Hier kann man nun das Ergebnis direkt ablesen:

x = -115/144   x = -13/8    x = 233/144
 0              1            2

Mit dem Programm TTMathe kannst du übrigens für jedes LGS diese Schritte verfolgen.

Folgendes Programm zeigt Dir, wie Du das zweidimensionale Feld koeff deklarierst und wie Du das LGS löst:
procedure TForm1.Button1Click(Sender: TObject);
  const n=3;
  var i0: integer;
      x: array[0..n-1] of extended;             //Die Lösungen x0, x1, x2, ...
      koeff: Array[0..n-1, 0..n] of extended; //Die Koeffizienten
  procedure VereinfacheRest(i: integer);
  var zeile, spalte: integer;
    d: extended;
  begin
    for zeile := 0 to n - 1 do if zeile <> i then Begin
        d := koeff[zeile, i] / koeff[i, i];
        if d <> 0 then Begin
          for spalte := 0 to n do if spalte <> i then
              koeff[zeile, spalte] := koeff[zeile, spalte] - d * koeff[i, spalte] else
              koeff[zeile, i] := 0; //=koeff[zeile,i]-koeff[zeile,i]/koeff[i,i]*koeff[i,i]
        End;
      End;
  end;
begin
  koeff[0,0] := 5;  koeff[0,1] := -7; koeff[0,2] := 1; koeff[0,3] := 9;
  koeff[1,0] := -3; koeff[1,1] := 2;  koeff[1,2] := 3; koeff[1,3] := 4;
  koeff[2,0] := 2;  koeff[2,1] := 3;  koeff[2,2] := 4; koeff[2,3] := 0;
  for i0 := 0 to n - 1 do VereinfacheRest(i0);
  for i0 := 0 to n - 1 do x[i0] := koeff[i0, n]/koeff[i0,i0] ;
  for i0 := 0 to n - 1 do
    showmessage('x' + intToStr(i0+1)+ ' = ' + FloatToStr(x[i0]));
end;
Beispiel 15.2 Der allgemeine Gauß'sche Algorithmus.
Diesmal verwenden wir ein dynamische Array. Das zweidimensionale Array wird deklariert durch den Typ:

type TarrayOfArrayOfExtended = array of array of extended;

Falls koeff[i,i] = 0 ist müssen noch die Spalten vertauscht werden. Das geschieht mit der Permutation p und Umkehrpermutation q.
Hier musst Du noch die Untit MathMohr miteinbinden. (oder die Entsprechenden Prozeduren durch eigene ersetzten.)
function IsInteger(const x: extended; eps: extended): boolean;
begin //eps globale Variable. Zum Beispiel eps = 1E-9
  result := frac(abs(x) + eps) < eps*2;
end;

FUNCTION ggtInt(a,b:longint):longint;
  begin if b=0 then result:=a
        else result:=ggtInt(b,a mod b);
  end;

FUNCTION ggtReal(a,b:Extended):Extended;
  begin if (a < maxlongint) and (b < maxlongint) then
           result:=ggTInt(round(a),round(b))
        else Begin
          if abs(b) < 0.5 then result:=a
          else result:=ggtReal(b,a-b*int(a/b));
        end;
end;

FUNCTION kgVReal(a, b: extended): extended;
begin result := a * b / ggTReal(a, b) end;

procedure LSG(n: integer; var aa: TarrayOfArrayOfExtended; var xx: array of extended);
             //Arrray aa[0..n,0..n+1]
var i0, j0,j0max: integer; //n Unbekannte bzw. Gleichungen
  p, q: array of integer; //permutation q=inv_p
  procedure invers; //erzeugt q = inverse Permutation von p
  var u, v: integer;
  begin
    for u := 0 to n - 1 do
      for v := 0 to n - 1 do
        if p[u] = v then q[v] := u;
  end;
  procedure tausche_sp(i, j: integer); // Spalten werden ausgetauscht
  var u, k: integer;
    x: extended;
  begin
    for u := 0 to n - 1 do Begin
      x := aa[u, i];
      aa[u, i] := aa[u, j];
      aa[u, j] := x; //=altes aa[u,i]
    End;
    k := p[i];
    p[i] := p[j];
    p[j] := k; //altes p[i]
    invers;
  end;

  procedure macheZeileGanzzahlig(zeile: integer); //bis auf rechte Seite  aa[z,n+1]
  var k : integer;
    d, zae, nen: extended; //wird das kgV des Nenners
  begin
    try
      d := 1;
      for k := 0 to n - 1 do Begin
        if not ErmittleBruch(abs(aa[zeile, k]), zae, nen,g_eps) then exit;
        d := kgVReal(nen, d);
      End;
      for k := 0 to n do aa[zeile, k] := d * aa[zeile, k];
      //Jetzt noch kürzen
      if not isInteger(aa[zeile, 1], g_eps) then exit;
      d := round(aa[zeile, 1]);
      for k := 0 to n - 1 do Begin
        if not isInteger(aa[zeile, k], g_eps) then exit;
        if d = 0 then d := round(aa[zeile, k]); //falls aa[zeile,1..]=0
        if abs(aa[zeile, k]) > 0 then
          d := ggTReal(round(aa[zeile, k]), d);
      End;
      if d <> 0 then for k := 0 to n do aa[zeile, k] := aa[zeile, k] / d;
    except {dann halt nicht} end;
    if aa[zeile, zeile] < 0 then for k := 0 to n do aa[zeile, k] := -aa[zeile, k];
  end;

  procedure VereinfacheRest(i: integer);
  var zeile, spalte: integer;
    d: extended;
  begin
    for zeile := 0 to n - 1 do if zeile <> i then Begin
        d := aa[zeile, i] / aa[i, i];
        if d <> 0 then Begin
          for spalte := 0 to n do if spalte <> i then
              aa[zeile, spalte] := aa[zeile, spalte] - d * aa[i, spalte] else
              aa[zeile, i] := 0; //=aa[zeile,i]-aa[zeile,i]/aa[i,i]*aa[i,i]
          macheZeileGanzzahlig(zeile);
        End;
      End;
  end;

begin //Hauptprogramm
  setlength(p,n+1);
  setlength(q,n);
  for j0 := 0 to n - 1 do Begin
    p[j0] := j0;
    q[j0] := j0;
  End;
  for i0 := 0 to n - 1 do Begin
    j0max := i0;
    for j0 := i0 + 1 to n - 1 do if abs(aa[i0, j0]) > abs(aa[i0, j0max]) then j0max := j0;
    if aa[i0,j0max] = 0 then Begin
      showmessage('Keine eindeutige Lösng!');
      exit;
    End;
    VereinfacheRest(i0);
  end;
  for i0 := 0 to n - 1 do xx[p[i0]] := aa[i0, n] ;
end;
Die Funktione ermittle Bruch findet man in der unit jmbruch
procedure ArrayInMemo(n: integer; aa: TarrayOfArrayOfExtended; m:Tmemo);
  var i, j: integer;   //aa n zeilen und n + 1 Spalten
      s: string;
begin
  m.Lines.Clear;
  for i := 0 to n-1 do Begin
    s := '';
    for j := 0 to n do
      s := s + ' ' + floatToStr(aa[i,j]);
    m.Lines.Add(s);
  End;
end;

procedure TForm1.Button1Click(Sender: TObject);
  var n, i: integer;
      koeff: TarrayOfArrayOfExtended;
      a: array of extended;
begin
  n := 3; // von 0 bis 2
  setlength(a, n);
  setlength(koeff, length(a),length(a) + 1);
  koeff[0,0] := 5;  koeff[0,1] := -7; koeff[0,2] := 1; koeff[0,3] := 9;
  koeff[1,0] := -3; koeff[1,1] := 2;  koeff[1,2] := 3; koeff[1,3] := 4;
  koeff[2,0] := 2;  koeff[2,1] := 3;  koeff[2,2] := 4; koeff[2,3] := 0;
  ArrayInMemo(length(a),Koeff,memo1);
  LSG(length(a),koeff,a);
  memo1.lines.add('Lösung');
  for i := 0 to n - 1 do
    memo1.lines.add(ReellZuBruch(a[i]));
end;
Hinweis: Downloadseite "LGS"

Aufgabe 15.1

Siehe dazu auch: Vollständige Induktion und Errechnung der Summenformel nach der Faulhaberformel.

Schreibe ein Programm, das mit Hilfe der Funktion LGS von Beispiel 12.2 eine geschlossene Formel für die folgenden Summen ermittelt:

a) S(n,k) = 12 + 22 + 32 + ... + n2
b) Allgemein für k = 0, 1, 2, 3, ... S(n,k) = 1k + 2k + 3k + ... + nk

Hinweis zu a) Ansatz:

               0      1      2      3
  S(n,2) = a0*n + a1*n + a2*n + a3*n



  Punktprobe mit n=0, 1 ,2 und 3 ergibt folgendes LGS:


               0      1      2      3
  S(0,2) = a0*0 + a1*0 + a2*0 + a3*0


               0      1      2      3
  S(1,2) = a0*1 + a1*1 + a2*1 + a3*1


               0      1      2      3
  S(2,2) = a0*2 + a1*2 + a2*2 + a3*2


               0      1      2      3
  S(3,2) = a0*3 + a1*3 + a2*3 + a3*3



Hinweis zu b) Ansatz


         S(n,k) = 1 + 2^k + 3^k + 4^k + ... + n^k

              0       1       2       2               k+1
         = a*n + a * n + a * n + a * n + ... + a   * n        (k+2 Unbekannte)
            0     1       2       3             k+1


              0      1      2             k+1                  0             1
  S(0,k) = a*0 + a *0 + a *0 + ... +a   *0       koeff[0,0] = 0  koeff[0,1]=0  ...
            0     1      2           k+1


              0     1     2             k+1                   0              1
  S(1,k) = a*1 + a*1 + a*1  + ... + a* 1        koeff[1,0] = 1 koeff[1,1] = 1  ...
            0     1     2            k+1


              0             2              k+1                0              1
 S(2,k) = a*2  + a *2 + a *2  + ... + a   2    koeff[2,0] = 2 koeff[2,1] = 2  ...
            0      1      2             k+1


    ...

                     0          1         2
  S(k+1,k) = a *(k+1) + a *(k+1) + a (k+1)  + ... + a     (k+1)^(k+1)
              0          1          2                 k + 1


                                0                     1
            koeff[k+1,0] = (k+1)  koeff[k+1,1] = (k+1)  ...



  Also: koeff[i,j] = i^j  (i=0 ... k+1; j=0 ... k+1; k+1 Gleichungen mit k+1 Unbekannten


Lösung