Vorlesung 9: Anwendungen aus der linearen Algebra

In der numerischen linearen Algebra verwendet man sowohl direkten Methoden als auch iterative Methoden . In dieser Vorlesung beschaeftigen wir uns mit direkten Verfahren.

Called Functions

bslashtx
lutx

Loesung linearer Gleichunsgssysteme

Wir betrachten zunaechst die Loesung linearer Gleichungssysteme unter Verwendung des Gaussschen Eliminationsverfahrens.

Betrachte Ax=b mit einer quadratischen nxn Matrix A, det A ungleich 0.

Beispiel 1:

10x - 7y = 7

-3x + 2y +6z = 4

5x - y + 5z = 6

A = [10, -7, 0;
    -3, 2, 6;
    5, -1, 5];
b = [7;4;6];
x = A\b

Beispiel 2: Loese A x_i = b_i fuer verschiedene rechte Seiten b_1,..., b_k

A = magic(5);
B =  ones(5,3);
B(2,2) = 0;
B(4,3) = 0;
X = A\B

X =

    0.0154   -0.0358    0.0142
    0.0154    0.0527    0.0027
    0.0154    0.0123    0.0123
    0.0154    0.0219   -0.0281
    0.0154    0.0104    0.0604

Probe:

A*X

ans =

    1.0000    1.0000    1.0000
    1.0000   -0.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    0.0000
    1.0000    1.0000    1.0000

Oftmals sind die Vektoren b_i der rechten Seite nicht bereits zu Beginn der Rechnung bekannt, sondern werden erst im Verlauf der Rechnung bestimmt. In solch einer Situation ist die Verwendung von A\ ineffizient, da viele Rechenschritte mehrfach wiederholt werden. Stattsessen sollte man zunaechst eine LU-Zerlegung der Matrix A bestimmen und mit dieser Zerlegung die verschiedenen Gleichungssysteme loesen. Wir betrachten noch einmal das obige Beispiel in etwas abgeaenderter Form.

Naive Implementation

clear all;
A = magic(5);
B = ones(5,3);
B(2,2) = 0;
B(4,3) = 0;
X = zeros(5,3);

tic
for k=1:3
    b = B(:,k);
    X(:,k) = A\b;
end
toc

Elapsed time is 0.000130 seconds.

Implementation unter Verwendung der LU-Zerlegung

clear all;
A = magic(5);
B = ones(5,3);
B(2,2) = 0;
B(4,3) = 0;
X = zeros(5,3);

tic
[L,U,P] = lu(A);
for k=1:3
    b = B(:,k);
    X(:,k) = U\(L\(P*b));
end
toc

Elapsed time is 0.000190 seconds.

Inverse Iteration

Algorithmus zur Bestimmung des Eigenvektors zum betragskleinsten Eigenwert. (Bei diesem Algorithmus loest man eine Folge von linearen Gleichungssystemen mit gleicher Matrix A und unterschiedlichen rechten Seiten b_i. Die Vektoren b_i werden erst im Verlauf der Rechnung bestimmt.)

Naive Implementation

clear all;
rng('default');
tic
A = randn(400)/20;
v = randn(400,1);
v = v/norm(v); % Anfangsnaeherung
for k=1:50
    v = A\v;
    v = v/norm(v);
end
lambda = (v'*A*v) / (v'*v); % Berechne zugehoerigen Eigenwert
toc
norm(A*v-lambda*v)  % Probe

Elapsed time is 0.110221 seconds.

ans =

   7.3140e-14

Implementation unter Verwendung der LU-Zerlegung: PA = LU

clear all;
rng('default');
tic
A = randn(400)/20;
[L,U,P] = lu(A);
v = randn(400,1);
v = v/norm(v); % Anfangsnaeherung
for k=1:50
    v = U \ (L\(P*v));
    v = v/norm(v);
end
lambda = (v'*A*v) / (v'*v); % Berechne zugehoerigen Eigenwert
toc
norm(A*v-lambda*v)  % Probe

Elapsed time is 0.015812 seconds.

ans =

   7.3140e-14

Permutationen

Eine Permutationsmatrix P ist eine quadratische Matrix, die in jeder Zeile und jeder Spalte genau eine 1 hat, alle anderen Elemente sind 0. Linksmultiplikation mit einer Permutationsmatrix fuehrt zu einem Zeilentausch, Rechtsmultiplikation mit einer Permutationsmatrix fuehrt zu einem Spaltentausch.

Beispiel:

A = [10, -7, 0; -3, 2, 6; 5, -1, 5];
P = [1, 0, 0; 0, 0, 1; 0, 1, 0];
A
P*A
A*P

A =

    10    -7     0
    -3     2     6
     5    -1     5


ans =

    10    -7     0
     5    -1     5
    -3     2     6


ans =

    10     0    -7
    -3     6     2
     5     5    -1

In Matlab koennen wir auch einen Permutationsvektor verwenden.

p = [1, 3, 2];

A(p,:) liefert das gleiche Ergebnis wir P*A.

A(p,:)

ans =

    10    -7     0
     5    -1     5
    -3     2     6

A(:,p) liefert das gleiche Ergebnis wie A*P.

A(:,p)

ans =

    10     0    -7
    -3     6     2
     5     5    -1

Matlab Implementation der LU Zerlegung

lutx ist eine Matlab Implementation der LU Zerlegung.

Autor: C.Moler

clear all
A = [10, -7, 0;
    -3, 2, 6;
    5, -1, 5];
[L,U,p]=lutx(A)

L =

    1.0000         0         0
    0.5000    1.0000         0
   -0.3000   -0.0400    1.0000


U =

   10.0000   -7.0000         0
         0    2.5000    5.0000
         0         0    6.2000


p =

     1
     3
     2

Matlab Implementation von backslash

bslashtx ist eine Matlab Implementation von backslash.

Autor: C. Moler

bslashtx prueft zunaechst ob es sich um eine untere Dreiecksmatrix, eine obere Dreiecksmatrix oder um eine symmetrisch positiv definite Matrix handelt. In diesen Spezialfaellen koennen effizientere Algorithmen verwendet werden.

Falls keiner dieser Spezialfaelle vorliegt, dann wird unter Verwendung von lutx eine LU Zerlegung mit Permutationsmatrix berechnet.

bslashtx enthaelt subfunctions zur Loesung von linearen Gleichungssystmen mit unterer oder oberer Dreiecksmatrix.

clear all;
A = [10, -7, 0;
    -3, 2, 6;
    5, -1, 5];
b = [7;4;6];
x = bslashtx(A,b)

function x = bslashtx(A,b)
% BSLASHTX  Solve linear system (backslash)
% x = bslashtx(A,b) solves A*x = b

[n,n] = size(A);
if isequal(triu(A,1),zeros(n,n))
   % Lower triangular
   x = forward(A,b);
   return
elseif isequal(tril(A,-1),zeros(n,n))
   % Upper triangular
   x = backsubs(A,b);
   return
elseif isequal(A,A')
   [R,fail] = chol(A);
   if ~fail
      % Positive definite
      y = forward(R',b);
      x = backsubs(R,y);
      return
   end
end

% Triangular factorization
[L,U,p] = lutx(A);

% Permutation and forward elimination
y = forward(L,b(p));

% Back substitution
x = backsubs(U,y);


% ------------------------------

function x = forward(L,x)
% FORWARD. Forward elimination.
% For lower triangular L, x = forward(L,b) solves L*x = b.
[n,n] = size(L);
x(1) = x(1)/L(1,1);
for k = 2:n
   j = 1:k-1;
   x(k) = (x(k) - L(k,j)*x(j))/L(k,k);
end


% ------------------------------

function x = backsubs(U,x)
% BACKSUBS.  Back substitution.
% For upper triangular U, x = backsubs(U,b) solves U*x = b.
[n,n] = size(U);
x(n) = x(n)/U(n,n);
for k = n-1:-1:1
   j = k+1:n;
   x(k) = (x(k) - U(k,j)*x(j))/U(k,k);
end

function [L,U,p] = lutx(A)
%LUTX  Triangular factorization, textbook version
%   [L,U,p] = lutx(A) produces a unit lower triangular matrix L,
%   an upper triangular matrix U, and a permutation vector p,
%   so that L*U = A(p,:)

[n,n] = size(A);
p = (1:n)';

for k = 1:n-1

   % Find index of largest element below diagonal in k-th column
   [r,m] = max(abs(A(k:n,k)));
   m = m+k-1;

   % Skip elimination if column is zero
   if (A(m,k) ~= 0)

      % Swap pivot row
      if (m ~= k)
         A([k m],:) = A([m k],:);
         p([k m]) = p([m k]);
      end

      % Compute multipliers
      i = k+1:n;
      A(i,k) = A(i,k)/A(k,k);

      % Update the remainder of the matrix
      j = k+1:n;
      A(i,j) = A(i,j) - A(i,k)*A(k,j);
   end
end

% Separate result
L = tril(A,-1) + eye(n,n);
U = triu(A);