function Aufg_6_3

close all
dbstop if error

%impl. Mittelpunkregel:
A1 = 1/2;
b1 = 1; c1=1/2;

%Crank-Nicolson Methode:
A2 = zeros(2,2);    A2(2,1)=1/2;    A2(2,2)=1/2;
b2 = [1/2;1/2]; c2=[0;1];

% Startwerte:
y0=[0.8;0];
t0=0;

% function handles
f = @(t,y) [y(2); -3*y(1)]; % rechte Seite der DGL
fy = @(t,y) [0 1; -3 0]; % f_y -- Ableitung der rechten Seite
y = @(t) [0.8*cos(sqrt(3)*t);-0.8*sqrt(3)*sin(sqrt(3)*t)]; % exakte Lösung des AWP im b-Teil

% Deklarationen:
f_err=zeros(10,2);
T=1;
y_von_1= y(1);

h = (1/2).^(1:10); % Schrittweiten

% Berechnung der Approximationen an y(1) mit den verschiednen
% Schrittweiten h(i), i=1,...,10:
for i = 1:10 
        
    % Berechnung der Approx. mit den verschiedenen Verfahren:
    [~,y1] = impl_Runge_Kutta_Verfahren(t0,y0,h(i),f,fy,T,A1,b1,c1);
    [~,y2] = impl_Runge_Kutta_Verfahren(t0,y0,h(i),f,fy,T,A2,b2,c2);
    
    %hier werden die relativen Fehler berechnet:
    f_err(i,1)= norm(y1(:,end)-y_von_1);
    f_err(i,2)= norm(y2(:,end)-y_von_1);
   
end

%Plot der Norm des absoluten Fehlers gegen die Schrittweite h (in loglog):
loglog(h,f_err(:,1),'-+',h,f_err(:,2),'-o',h,h.^2,':');
legend('impl. Mittelpunktregel','Crank-Nicolson','h^2','location','NorthWest');
title('Fehlerplot')
ylabel('Norm des absoluten Fehlers')
xlabel('Schrittweite h')

end


function [ti,yi] = impl_Runge_Kutta_Verfahren(t0,y0,h,f,fy,T,A,b,c)
% [ti,yi] diskreter Loesungsgraph des AWPs y=f(y), y(to)=y0
% auf dem Intervall [t0,T]
%
% A,b,c <--- Buchter-Parameter des (impl.) RKVs


% Initialisierung:
ti=t0:h:T; % Zeitgitter
s=max(size(b)); % s ist die Stufenzahl des Verfahrens
N=length(ti); % Anzahl der Zeitschritte
d=max(size(y0)); % Dimension des AWPs
yi=zeros(d,N); % Loesungsvektor (Approximationen an die Lsg. y)
F=zeros(d*s,1); % Vektor aus dem Newton-Verfahren vgl. Aufg. 5.1
AFy=zeros(d*s,d*s); % ein Ableitungsterm aus dem Newton-Verfahren vgl. A 5.1

yi(:,1)=y0; % Startwert

for i=1:N-1 % Iterationsschritte des impliziten RKVs
    
    Yi=kron(ones(s,1),yi(:,i)); % ein Vektor der nach Aufg. 5.1 in dem
                                % Newton-Verfahren auftaucht 
                                
    ETA=Yi; % Startwert des Newton Verfahrens
    
    continue_loop = 1; % Newton Schleife initialisieren
    count = 0;
    
    while continue_loop % Newton-Verfahren zur Bestimmung von ETA bzw. eta 
        
        % Auswerten von f in eta_1 .... eta_s und (A kron E_dxd) *
        % Auswertungen von fy bestimmen.
        for l = 1:s
            F((l-1)*d+1:l*d,1)=f(ti(i)+c(l)*h,ETA((l-1)*d+1:l*d));
            AFy(:,(l-1)*d+1:l*d) = kron(A(:,l),fy(ti(i)+c(l)*h,ETA((l-1)*d+1:l*d)));
        end
        % N' und N von Aufgabe 5.1 bestimmen
        Nprime = speye(d*s) - h*AFy;
        N = ETA - Yi - h*kron(A,speye(d))*F;
        % Newtonschritt amchen
        ETA = ETA - Nprime \ N;
        
        count = count +1;
                
        if norm(N)<10^(-8) % eigentlicher Newtonabbruch
            continue_loop = 0;
        end
        
        if count > 100 % alternativer Newtonabbruch
            error('Newton iteration seems to not converge!')
        end
        
    end
    
    eta=reshape(ETA,d,s); % ETA wird in eine Matrix eta mit den eta_j's in den Spalten umgeformt
        
    % Endauswertung eines RKV-Schrittes
    sum=0;
    for k=1:s
        sum=sum+b(k)*f(ti(i)+c(k)*h,eta(:,k));
    end
    yi(:,i+1)=yi(:,i)+h*sum;
    
end

end
