Musterlösung Aufgabe 24

Aufgabe 24: (LR-Faktorisierung von Bandmatrizen)

Das n × n-Gleichungssystem ( \begin{matrix}2 & −1 & & 0 \\ −1 & ... & ... \\ & ... & ... & −1 \\ 0 & & −1 & 2 \end{matrix} ) ( x_{1} ; ⋮ ; ⋮ ; x_{n} ) = ( b_{1} ; ⋮ ; ⋮ ; b_{n} ) ist zu lösen.

a) Finde die (exakten) Daten der LR-Zerlegung A = LR.

Anleitung: die die Bänder von L und R definierende Rekursion ist leicht lösbar.

b) Berechne (exakt) den Lösungsvektor x⃗ für b⃗ = (1, …, 1)^{T} .

Anleitung: die Rekursionen der Rücksubstitutionen sind für dieses b⃗ leicht lösbar.

c) Berechne numerisch die Komponente x_{10 000} für b⃗ = (sin(1), sin(2), …, sin(n))^{T}, n = 500 000. Welcher Speicherplatzbedarf fällt an? Notiere die Laufzeit, Quelltext des Programmes mitabgeben!

Achtung! Das Gleichungssystem ist schlecht konditioniert, relative Rundungsfehler werden um Faktoren der Größenordnung n^{2} verstärkt!

Musterlösung: a) Ansatz

LR = ( \begin{matrix}1 \\ l_{1} & ... \\ & ... & ... \\ & & l_{n−1} & 1 \end{matrix} ) ( \begin{matrix}r_{1} & −1 \\ & ... & ... \\ & & ... & −1 \\ & & & r_{n} \end{matrix} )

= ( \begin{matrix}r_{1} & −1 \\ l_{1} r_{1} & −l_{1} + r_{2} & −1 \\ & l_{2} r_{2} & −l_{2} + r_{3} & −1 \\ & & ... & ... & ... \\ & & & ... & ... & −1 \\ & & & & l_{n−1} r_{n−1} & −l_{n−1} + r_{n} \end{matrix} ) .

Durch Vergleich mit A folgt

r_{1} = 2 , l_{i} r_{i} = −1 (i = 1, …, n − 1) , −l_{i} + r_{i+1} = 2 (i = 2, …, n − 1) ,

also l_{i} = −1/r_{i} und

r_{i+1} = 2 + l_{i} = 2 − \frac{1}{r_{i}} , r_{1} = 2 .

Diese Rekursion liefert

r_{1} = 2 , r_{2} = 2 − \frac{1}{2} = \frac{3}{2} , r_{3} = 2 − \frac{2}{3} = \frac{4}{3} , ….

In der Tat ist

r_{i} = \frac{i + 1}{i} , i = 1, …, n , l_{i} = −\frac{1}{r_{i}} = −\frac{i}{i + 1}

die Lösung der Rekursion:

r_{i+1} = \frac{i + 2}{i + 1} =^{!} 2 − \frac{1}{r_{i}} = 2 − \frac{i}{i + 1} = \frac{i + 2}{i + 1} .

b) Rücksubstitutionen LRx⃗ = b⃗. Zunächst Ly⃗ = b⃗:

( \begin{matrix}1 \\ l_{1} & ... \\ & ... & ... \\ & & l_{n−1} & 1 \end{matrix} ) ( y_{1} ; ⋮ ; ⋮ ; y_{n} ) = ( b_{1} ; ⋮ ; ⋮ ; b_{n} )

⇒ y_{1} = b_{1} , y_{i} = b_{i} − l_{i−1} y_{i−1} = b_{i} + \frac{i − 1}{i} y_{i−1} (i = 2, …, n) .  (1)

Mit b_{i} = 1 ergeben sich aus

y_{1} = 1 , y_{i} = 1 + \frac{i − 1}{i} y_{i−1} (i = 2, …, n)

die Werte

y_{1} = 1 , y_{2} = 1 + \frac{1}{2} = \frac{3}{2} , y_{3} = 1 + \frac{2}{3} \frac{3}{2} = \frac{4}{2} , y_{4} = 1 + \frac{3}{4} \frac{4}{2} = \frac{5}{2} , · · · .

Die offensichtliche Lösung der Rekursion ist y_{i} = (i + 1)/2 :

y_{i} = \frac{i + 1}{2} =^{!} 1 + \frac{i − 1}{i} y_{i−1} = 1 + \frac{i − 1}{i} \frac{i}{2} = \frac{i + 1}{2} .

Zweite Rücksubstitution Rx⃗ = y⃗:

( \begin{matrix}r_{1} & −1 \\ & ... & ... \\ & & ... & −1 \\ & & & r_{n} \end{matrix} ) ( x_{1} ; ⋮ ; ⋮ ; x_{n} ) = ( y_{1} ; ⋮ ; ⋮ ; y_{n} )

liefert

x_{n} = \frac{y_{n}}{r_{n}} = \frac{n}{n + 1} y_{n} , x_{i} = \frac{y_{i} + x_{i+1}}{r_{i}} = \frac{i}{i + 1} (y_{i} + x_{i+1}) (i = n − 1, …, 1) , (2)

also

x_{n} = \frac{n + 1}{2} \frac{n}{n + 1} = \frac{n}{2} , x_{i} = (\frac{i + 1}{2} + x_{i+1}) \frac{i}{i + 1} = \frac{i}{2} + \frac{i}{i + 1} x_{i+1} .

Für z_{i} := x_{i}/i folgt die einfachere Rekursion

z_{n} = \frac{1}{2} , z_{i} = \frac{1}{2} + z_{i+1} (i = n − 1, …, 1)

mit der offensichtlichen Lösung

z_{i} = \frac{n + 1 − i}{2} .

Endergebnis:

x_{i} = i z_{i} = \frac{i (n + 1 − i)}{2} .

c) Für allgemeines b⃗ wird die Lösung mittels (1) und (2) durch

y_{1} := b_{1} , for i = 2..n do y_{i} := b_{i} + \frac{i − 1}{i} y_{i−1} ;

x_{n} := \frac{n}{n + 1} y_{n} , for i = n − 1..1 do x_{i} := \frac{i}{i + 1} (y_{i} + x_{i+1})

geliefert. Speichert man erst y⃗, dann x⃗ in b⃗, so liefert

for i = 2 .. n do b[i] := b[i] + (i − 1)/i ∗ b[i − 1] ;
b[n] := n/(n + 1) ∗ b[n] ;
for i = n − 1 .. 1 do b[i] := i/(i + 1) ∗ (b[i] + b[i + 1]) ;

die Lösung im Vektor b⃗. Der Speicherplatzbedarf besteht aus einem Vektor b⃗. Bei 64-bit(=8-Byte)-floats (16 Dezimalstellen) mit n = 500 000 sind dies 4 MByte.

MuPAD-Lösung mit verschiedenen Werten von n und verschiedenen Maschinengenauigkeiten:

   *----*    MuPAD 1.4.0 --- Multi Processing Algebra Data Tool
  /|   /|
 *----* |    Copyright (c) 1992-97 by B. Fuchssteiner, Automath
 | *--|-*    University of Paderborn. All rights reserved.
 |/   |/
 *----*      ----------- Developers NSB Version -----25.11.97--

>> for n in [5000,50000,500000] do
&>  for digits in [7,16,32] do
&>    DIGITS:=digits;
&>    st:=time():
&>    b:=array(1..n,[ float(sin(i)) $ hold(i)=1..n] ):
&>    for i from 2 to n do b[i]:=b[i]+(i-1)/i*b[i-1]: end_for:
&>    b[n]:=n/(n+1)*b[n];
&>    for i from n-1 downto 1 do b[i]:=i/(i+1)*(b[i]+b[i+1]); end_for:
&>    DIGITS:=10;
&>    print("n=".n,"DIGITS=".digits,"Ergebnis", b[n/50],"Zeit", (time()-st)*sec/1000.0);
&>  end_for;
&> end_for:
 "n=5000",   "DIGITS= 7", "Ergebnis", -0.5418536113, "Zeit",    7.56sec
 "n=5000",   "DIGITS=16", "Ergebnis", -0.5419803388, "Zeit",    8.29sec
 "n=5000",   "DIGITS=32", "Ergebnis", -0.5419803388, "Zeit",   10.54sec
 "n=50000",  "DIGITS= 7", "Ergebnis",  0.8984039626, "Zeit",   75.37sec
 "n=50000",  "DIGITS=16", "Ergebnis",  0.9114517129, "Zeit",   84.07sec
 "n=50000",  "DIGITS=32", "Ergebnis",  0.9114517129, "Zeit",  106.51sec
 "n=500000", "DIGITS= 7", "Ergebnis",  0.08038886289,"Zeit",  756.83sec
 "n=500000", "DIGITS=16", "Ergebnis", -0.3164850488, "Zeit",  829.65sec
 "n=500000", "DIGITS=32", "Ergebnis", -0.3164850491, "Zeit", 1051.66sec

C-Lösung mit double-Variablen (16-Stellen-Arithmetik):

#include <stdio.h>
#include <math.h>
#define  n 500000
main()
{ double b[n+1];
  int    i;
  for (i=1; i<=n; i++)   b[i]= sin( (double) i);
  for (i=2; i<=n; i++)   b[i]= b[i]+((i-1)*b[i-1])/i;
  b[n] = (n*b[n])/(n+1);
  for (i=n-1; i>=1; i--) b[i]= (i*(b[i]+b[i+1]))/(i+1);
  printf("%+.10f\n",b[10000]);
}
Ergebnis: -0.3164850492     Zeit:  1.46u 0.06s 0:01.55 98.0%