H06 – Iterative Lösung (Gradientenverfahren, CG)

 

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. Schreiben Sie ein m-File („u61.m“), in dem Sie das Gleichungssystem mit dem Gradientenverfahren und dem CG-Verfahren (jeweils ohne Vorkonditionierung, d.h. B = I) lösen. Verwenden Sie einen zufälligen Startvektor {x_0}.
Führen Sie die Rechnungen für N = 10,20,40,80 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? Wie erklären Sie sich das Ergebnis?

Lösung

u61.m

clear;
tic

versuche = 10;
N = [10 20 40 80];
IT = zeros (10,2);

for n = 1 : length(N)
    for versuch = 1 : versuche
        v = -1 * ones(1, N(n) - 1)';
        v1 = 2 * ones(1, N(n) - 1)';
        A = spdiags([v, v1, v], -1 : 1, N(n) - 1, N(n) - 1);
        A = full(A); % mit der vollständigen Matrix laufen die Verfahren schneller
        b = ones(N(n) - 1, 1);
        B = eye(N(n) - 1);

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

        % Gradientenverfahren

        xk = x0;
        r0 = b - A * x0;
        rk = r0;
        s0 = B * r0;
        sk = s0;
        v0 = A * s0;
        vk = v0;

        if r0 == 0
            disp('x0 ist exakt')
            return
        end

        k = 0;
        while norm(b - A * xk) >= exit
            ak = (rk' * sk) / (vk' * sk);
            xk = xk + ak * sk;
            rk = rk - ak * vk;
            sk = B * rk;
            vk = A * sk;
            k = k + 1;
        end

        % CG-Verfahren

        xk_cg = x0;
        r0_cg = b - A * x0;
        rk_cg = r0_cg;
        s0_cg = B * r0_cg;
        sk_cg = s0_cg;
        w0_cg = s0_cg;
        y0_cg = r0_cg' * w0_cg;
        yk_cg = y0_cg;
        k_cg = 0;

        while norm(b - A * xk_cg) >= exit
            vk_cg = A * sk_cg;
            dk_cg = sk_cg' * vk_cg;
            ak_cg = yk_cg/dk_cg;
            xk_cg = xk_cg + ak_cg * sk_cg;
            rk_cg = rk_cg - ak_cg * vk_cg;
            wk_cg = B * rk_cg;
            yk2_cg = rk_cg' * wk_cg;
            betak_cg = yk2_cg/yk_cg;
            sk_cg = wk_cg + betak_cg * sk_cg;
            yk_cg = yk2_cg;
            k_cg = k_cg + 1;
        end

        IT(versuch,1) = k;
        IT(versuch,2) = k_cg;
    end
    disp(['Iterationsschritte Grad.: ',num2str(sum(IT(:,1)/length(IT(:,1))))])
    disp(['Iterationsschritte CG: ',num2str(sum(IT(:,2)/length(IT(:,2))))])
end 

toc % Rechenzeit: 2.669970 Sekunden

% Iterationsschritte:
%
% N     Grad	CG
% 10    122.8   9
% 20    507.3   18.5
% 40    2026.4  33.1
% 80    7973.2  55.7
%
% Grad: quadratisch
% CG:   linear -> schneller, weil es nicht in Richtung des Gradienten geht,
%       sondern eine verbesserte Suchrichtung hat. Bessere Kondition