Joachim Mohr Mathematik Musik Delphi
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