H05 – Lösen eines Gleichungssystems (Jacobi, Gauß-Seidel, SOR)

 

Es ist das Gleichungssystem Ax = b mit

A = \left( {\begin{array}{*{20}{c}} 2 & {-1} & {} & {} \\ {-1} & 2 & {-1} & {} \\ {} & \ddots & \ddots & \ddots \\ {} & {} & {-1} & 2 \\ \end{array} } \right) \in {\mathbb{R}^{\left( {N-1} \right) \times \left( {N-1} \right)}},\quad b = \left( {\begin{array}{*{20}{c}} 1 \\ \vdots \\ 1 \\ \end{array} } \right) \in {\mathbb{R}^{\left( {N-1} \right)}}

zu lösen. Die Matrix A hat die Eigenwerte

{\lambda _i} = 2-2\cos \left( {\frac{i}{N}\pi } \right),\quad i = 1, \ldots ,N-1.

  1. Bestimmen Sie den Spektralradius der Iterationsmatrix {M_J} für das Jacobi-Verfahren
  2. Für den optimalen Relaxationsparameter {\omega _{opt}} des SOR-Verfahrens gilt

    {\omega _{opt}} = \frac{2}{{1+\sqrt {1-\rho {{\left( {{M_J}} \right)}^2}} }}.

    Bestimmen Sie {\omega _{opt}}.

  3. Schreiben Sie ein m-File („u51.m“), in dem Sie das Gleichungssystem mit dem Jacobi-, Gauß-Seidel- und SOR-Verfahren lösen. Verwenden Sie einen zufälligen Startvektor {x_0}. Führen Sie das SOR-Verfahren dabei einmal mit dem in b) bestimmten optimalen Relaxationsparameter aus sowie mit einem beliebigen weiteren Wert \omega \in \left[ {0,2} \right]. Führen Sie die Rechnungen für N = 50 jeweils zehnmal durch und bestimmen Sie für jedes der Verfahren die durchschnittliche Anzahl an Iterationen, die notwendig ist, um die Bedingung

    {\left\| {b-Ax} \right\|_2} < {10^{-3}} \cdot {\left\| {b-A{x_0}} \right\|_2}

    zu erfüllen. Schreiben Sie das Ergebnis als Kommentar in Ihr m-File. Was fällt Ihnen bei den Iterationszahlen auf?

Hinweise:

In der Lösung zu Aufgabe c) darf der Backslash-Operator (“\”) nicht vorkommen. Orientieren Sie sich stattdessen an der komponentenweisen Darstellung der Verfahren.

Lösung

a )

Für die Iterationsmatrix des Jacobi-Verfahrens gilt:

{M_J} = I-{D^{-1}}A

= \left( {\begin{array}{*{20}{c}} 1 & {} & {} \\ {} & \ddots & {} \\ {} & {} & 1 \\ \end{array} } \right)-\left( {\begin{array}{*{20}{c}} {\frac{1}{2}} & {} & {} \\ {} & \ddots & {} \\ {} & {} & {\frac{1}{2}} \\ \end{array} } \right)\left( {\begin{array}{*{20}{c}} 2 & {-1} & {} & {} \\ {-1} & 2 & {-1} & {} \\ {} & \ddots & \ddots & \ddots \\ {} & {} & {-1} & 2 \\ \end{array} } \right)

= \left( {\begin{array}{*{20}{c}} 1 & {} & {} \\ {} & \ddots & {} \\ {} & {} & 1 \\ \end{array} } \right)-\left( {\begin{array}{*{20}{c}} 1 & {-\frac{1}{2}} & {} & {} \\ {-\frac{1}{2}} & 1 & {-\frac{1}{2}} & {} \\ {} & \ddots & \ddots & \ddots \\ {} & {} & {-\frac{1}{2}} & 1 \\ \end{array} } \right) = \left( {\begin{array}{*{20}{c}} 0 & {-\frac{1}{2}} & 0 & {} & {} \\ {-\frac{1}{2}} & 0 & {-\frac{1}{2}} & 0 & {} \\ {} & \ddots & \ddots & \ddots & {} \\ {} & {} & {-\frac{1}{2}} & 0 & {-\frac{1}{2}} \\ {} & {} & {} & {-\frac{1}{2}} & 0 \\ \end{array} } \right)

Wir können die Iterationsmatrix auch schreiben als:

{M_J} = I-{D^{-1}}A = I-\frac{1}{2}A

Daher ergeben sich die Eigenwerte aus den Eigenwerten von A:

{\mu _i} = 1-\frac{1}{2}{\lambda _i} = 1-1+\cos \left( {\frac{i}{N}\pi } \right) = \cos \left( {\frac{i}{N}\pi } \right)

Dieser Wert ist für i = 1 am größten. Der Spektralradius ist also:

\rho \left( {{M_J}} \right) = \cos \left( {\frac{1}{N}\pi } \right)

b )

Wir können nun den Spektralradius berechnen:

{\omega _{opt}} = \frac{2}{{1+\sqrt {1-\rho {{\left( {{M_J}} \right)}^2}} }} = \frac{2}{{1+\sqrt {1-{{\left( {\cos \left( {\frac{1}{N}\pi } \right)} \right)}^2}} }}

Für die 49×49 Matrix erhalten wir:

{\omega _{opt}} = \frac{2}{{1+\sqrt {1-{{\left( {\cos \left( {\frac{1}{{50}}\pi } \right)} \right)}^2}} }} = 1.88183838983223

c )

Matlab-Code:

N = 50; 		% Abmessung der quadratischen Matrix ist (N-1)
versuche = 10; 	% Anzahl der Versuche mit verschiedenen Zufallsvektoren

% Aufstellen des Gleichungssystems
v = -1 * ones(1, N - 1)';
v1 = 2 * ones(1, N - 1)';
A = spdiags([v, v1, v], -1 : 1, N - 1, N - 1);
A = full(A); % mit der vollständigen Matrix laufen die Verfahren schneller
b = ones(N - 1, 1);

% Berechnung der Iterationsmatrix (Jacobi) und des Spektralradius
v = -.5 * ones(1, N - 1)';
Mj = spdiags([v, v], -1 : 2: 1, N - 1, N - 1);
e = eig(Mj);
r = max(abs(e));

% Berechnung des optimalen Relaxationsparameters (SOR)
wOpt = 2 / (1 + sqrt(1 - r * r));

% Vektoren für die Anzahl der Iterationsschritte
itJacobi = zeros(versuche, 1);
itGSeidel = zeros(versuche, 1);
itSorOpt = zeros(versuche, 1); % optimal
itSorNO = zeros(versuche, 1);  % nicht optimal

for versuch = 1 : versuche
    % Zufälliger Startvektor und Abbruchkriterium
    x0 = rand(N - 1, 1);
    exit = 0.001 * norm((b - A * x0));

    % Jacobi-Verfahren Komponentenweise
    x1 = x0;
    x2 = x0;
    while norm(b - A * x2) >= exit
        for i = 1 : N - 1
            x2(i) = 1/A(i,i) * (b(i) - sum(A(i, 1:i-1)' .* x1(1:i-1)) ...
                                     - sum(A(i, i+1:N-1)' .* x1(i+1:N-1)));
        end
        itJacobi(versuch) = itJacobi(versuch) + 1;
        x1 = x2;
    end % end Jacobi

    % Gauß-Seidel-Verfahren Komponentenweise
    x1 = x0;
    x2 = x0;
    while norm(b - A * x2) >= exit
        for i = 1 : N - 1
            x2(i) = 1/A(i,i) * (b(i) - sum(A(i, 1:i-1)' .* x2(1:i-1)) ...
                                     - sum(A(i, i+1:N-1)' .* x1(i+1:N-1)));
        end
        itGSeidel(versuch) = itGSeidel(versuch) + 1;
        x1 = x2;
    end % end Gauß-Seidel   

    % SOR-Verfahren Komponentenweise, w_opt = 1.88183838983223
    x1 = x0;
    x2 = x0;
    while norm(b - A * x2) >= exit
        for i = 1 : N - 1
            x2(i) = x1(i) + wOpt * 1/A(i,i) * (b(i) ...
                                  - sum(A(i, 1:i-1)' .* x2(1:i-1)) ...
                                  - sum(A(i, i:N-1)' .* x1(i:N-1)));
        end
        itSorOpt(versuch) = itSorOpt(versuch) + 1;
        x1 = x2;
    end % end SOR 

    % SOR-Verfahren Komponentenweise, w = 0.8
    x1 = x0;
    x2 = x0;
    w = 0.8;
    while norm(b - A * x2) >= exit
        for i = 1 : N - 1
            x2(i) = x1(i) + w * 1/A(i,i) * (b(i) ...
                                  - sum(A(i, 1:i-1)' .* x2(1:i-1)) ...
                                  - sum(A(i, i:N-1)' .* x1(i:N-1)));
        end
        itSorNO(versuch) = itSorNO(versuch) + 1;
        x1 = x2;
    end % end SOR
end % end versuch

itJacobiMed = sum(itJacobi) / versuche    % 3359.6
itGSeidelMed = sum(itGSeidel) / versuche  % 1678.3
itSorOptMed = sum(itSorOpt) / versuche    % 84.7   bei w = w_opt
itSorNOMed = sum(itSorNO) / versuche      % 2517.5 bei w = 0.8

toc % 45.360618 Sekunden

% Was fällt mir bei den Iterationszahlen auf?
%
%  n    Jacobi   Gauß-Seidel  SOR opt   SOR2
%  5       32       16          8         25
% 10      133       67         17        101
% 15      304      153         25        230
% 20      539      267         34        401
%
% Jacobi: nötige Iterationsschritte steigen quadratisch
%
% Gauß-Seidel: Iterationsschritte vervierfachen sich bei Verdopplung von n
%
% SOR opt: Iterationsschritte steigen linear
%
% SOR2: Iterationsschritte vervierfachen sich bei Verdopplung von n