Es ist das Gleichungssystem mit

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. ) lösen. Verwenden Sie einen zufälligen Startvektor .
Führen Sie die Rechnungen für jeweils zehnmal durch und bestimmen Sie für jedes der Verfahren die durchschnittliche Anzahl an Iterationen, die notwendig ist, um die Bedingung

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