% Aufgabe 3
% Bastian Groß
% Universität Trier
% 27. Oktober 2011
clc; clear all;

%Parameter
y0 = 1;
T  = 10;

d  = 13;

f_exact = inline('exp(-t)+t','t');
f_exp   = inline('x - y + 1','x','y');
x = 0;

for m = 1 : d
    %Rechengrößen
    n  = 2^m;
    h  = T/n;
    xa = x;
    x  = 0:h:T;
    
    f_rk(1,m)         = y0; 
    f_richardson(1,m) = y0;
    
    for i = 1 : n
        k1 = f_exp(x(i),f_rk(i,m));
        k2 = f_exp(x(i)+h/2,f_rk(i,m)+k1*h/2);
        k3 = f_exp(x(i)+h,f_rk(i,m)-k1*h+2*h*k2);
        f_rk(i+1,m) = f_rk(i,m) + (h/6)*(k1 + 4*k2 + k3); %Runge Kutta Fehlberg Verfahren 3ter Ordnung
    end
    if m > 1
            for e = 2 : n/2 + 1
                f_richardson(e,m-1) = f_rk(2*e-1,m) + (f_rk(2*e-1,m)-f_rk(e,m-1))/(2^3-1); % Richardson Korrektur (Extrapolation)
            end
    end
    
    error_rk(m) = max(abs(f_rk(:,m) - f_exact(x)')); 
    dt(m)          = h;
    M(m)           = m;
    if m > 1
        error_rich(m-1)   = max(abs(f_richardson(:,m-1) - f_exact(xa)')); 
        dt2(m-1)          = h;
        M2(m-1)           = m;
    end

        
    if m==2
        hold on; title('Schrittweiten h = 1/4 und h = 1/8')
        plot(x,f_rk(:,m),'b');      
    elseif m==3
        plot(x,f_rk(:,m),'g'); 
        plot(xa,f_richardson(:,m-1),'r');
    end
end

plot(x,f_exact(x),'k'); legend('Runge Kutta 3. Ordnung 1/4','Runge Kutta 3. Ordnung 1/8','Richardson Korrektur 1/4','Exakte Lösung','Location','Best');

k_rk   = log(error_rk(d)/error_rk(d-1))/log(1/2);
k_rich = log(error_rich(d-1)/error_rich(d-2))/log(1/2);
 
fid = fopen('Programmieraufgabe_3.txt','w');
fprintf(fid,'\n');
fprintf(fid,'       Runge Kutta Verfahren   \n');
fprintf(fid,'\n');
fprintf(fid,'   h       error        error/h      error/h^2   error/h^3    error/h^4        \n');
fprintf(fid,'-----------------------------------------------------------------------\n');
fprintf(fid,' 1/2^%2.0f   %1.3e   %1.3e   %1.3e  %1.3e   %1.3e  \n',[M;error_rk;error_rk./dt;error_rk./(dt.^2);error_rk./(dt.^3);error_rk./(dt.^4)]);
fprintf(fid,'\n     Ordnung k = %1.3f     \n',k_rk);
fprintf(fid,'-----------------------------------------------------------------------\n');
fprintf(fid,'\n');
fprintf(fid,'       Richardson Korrektur   \n');
fprintf(fid,'\n');
fprintf(fid,'   h       error        error/h      error/h^2   error/h^3    error/h^4        \n');
fprintf(fid,'-----------------------------------------------------------------------\n');
fprintf(fid,' 1/2^%2.0f   %1.3e   %1.3e   %1.3e  %1.3e   %1.3e  \n',[M2;error_rich;error_rich./dt2;error_rich./(dt2.^2);error_rich./(dt2.^3);error_rich./(dt2.^4)]);
fprintf(fid,'\n     Ordnung k = %1.3f     \n',k_rich);
fprintf(fid,'-----------------------------------------------------------------------\n');
fprintf(fid,'\n');

