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

%Parameter
P0 = 2;
K  = 20;
a  = 0.1;
T  = 5;

d  = 13;

for m = 1 : d
    %Rechengrößen
    n  = 2^m;
    h  = 1/n;
    t  = 0:h:T;
    
    f_euler(1) = P0;
    f_heun(1)  = P0;
    f_heun3(1) = P0;
    
    
    for i = 1 : T*n
        f_euler(i+1) = f_euler(i) + h*f_ode(f_euler(i),a,K); %Expliziter Euler 
        f_heun(i+1)  = f_heun(i)  + 0.5*h*(f_ode(f_heun(i)+h*f_ode(f_heun(i),a,K),a,K)+f_ode(f_heun(i),a,K)); %Heun Verfahren
        f_heun3(i+1) = f_heun3(i) + 0.25*h*(f_ode(f_heun3(i),a,K) + 3*f_ode(f_heun3(i)+(2*h/3)*f_ode(f_heun3(i)+(h/3)*f_ode(f_heun3(i),a,K),a,K),a,K));%Heun Verfahren 3ter Ordnung
    end
    
    error_euler(m) = max(abs(f_euler - sprucebudworm_exact(t,a,K,P0))); %Approximationsfehler
    error_heun(m)  = max(abs(f_heun - sprucebudworm_exact(t,a,K,P0)));
    error_heun3(m) = max(abs(f_heun3 - sprucebudworm_exact(t,a,K,P0)));
    dt(m)          = h;
    M(m)           = m;
        
    if m==2
        subplot(2,1,1); hold on; title('Schrittweite h = 1/4')
        plot(t,f_euler,'b'); plot(t,f_heun,'g'); plot(t,f_heun3,'r');
    elseif m==10
       subplot(2,1,2); hold on; title('Schrittweite h = 1/1024')
       plot(t,f_euler,'b'); plot(t,f_heun,'g'); plot(t,f_heun3,'r');
    end
end

subplot(2,1,1); plot(t,sprucebudworm_exact(t,a,K,P0),'k');
legend('Euler','Heun','verb. Heun', 'Exaktes P(t)','Location','SouthEast');
subplot(2,1,2); plot(t,sprucebudworm_exact(t,a,K,P0),'k');
legend('Euler','Heun','verb. Heun', 'Exaktes P(t)','Location','SouthEast');

k_euler = log(error_euler(d)/error_euler(d-1))/log(1/2); % Konvergenzordnung
k_heun  = log(error_heun(d)/error_heun(d-1))/log(1/2);
k_heun3 = log(error_heun3(d)/error_heun3(d-1))/log(1/2);

fid = fopen('Programmieraufgabe_1.txt','w');
fprintf(fid,'\n');
fprintf(fid,'            spruce budworm   \n');
fprintf(fid,'\n');
fprintf(fid,'       Euler 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_euler;error_euler./dt;error_euler./(dt.^2);error_euler./(dt.^3);error_euler./(dt.^4)]);
fprintf(fid,'\n     Ordnung k = %1.3f     \n',k_euler);
fprintf(fid,'-----------------------------------------------------------------------\n');
fprintf(fid,'\n');
fprintf(fid,'       Heun 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_heun;error_heun./dt;error_heun./(dt.^2);error_heun./(dt.^3);error_heun./(dt.^4)]);
fprintf(fid,'\n     Ordnung k = %1.3f     \n',k_heun);
fprintf(fid,'-----------------------------------------------------------------------\n');
fprintf(fid,'\n');
fprintf(fid,'       verb. Heun 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_heun3;error_heun3./dt;error_heun3./(dt.^2);error_heun3./(dt.^3);error_heun3./(dt.^4)]);
fprintf(fid,'\n     Ordnung k = %1.3f     \n',k_heun3);
fclose(fid);

% hold off
% loglog(dt,error_euler,'b'); hold on; title('loglog Plot');
% loglog(dt,error_heun,'g'); loglog(dt,error_heun3,'r')
% legend('Euler','Heun','verb. Heun','Location','SouthEast');
% xlabel('Schrittweite'); ylabel('Fehler');