% Plottar felet i ENM 7.13 som funktion av h i en loglog-plott
% Joar Bagge 2016-11-11

function felplott
    % Gör en referenslösning med fint h
    href = 2^-15
    wref = solve(href);

    % Gör nu flera lösningen med olika h
    hvec = 2.^(-(0:15)) % 2^0, 2^-1, 2^-2, ..., 2^-15
    errvec = []; % vektor för felen
    for h=hvec
        w = solve(h);
        % Vi mäter felet i sista tidpunkten (därav 'end' nedan).
        % Vi tar hänsyn till alla 4 komponenter i lösningsvektorn (därav ':').
        err = norm(w(:,end)-wref(:,end));
        errvec = [errvec, err];
    end

    % Plotta resultaten i en loglog-plott
    clf
    loglog(hvec, errvec, '-*'), hold on
    loglog(hvec, 30*hvec, 'k--') % linjär referenslinje e ~ h
    xlabel('h')
    ylabel('Fel')
    % Vi förväntar oss att felkurvan ska vara en rät linje med samma
    % lutning som referenslinjen. Detta stämmer inte för stora h eftersom
    % resultaten för noggrannhetsordningen bara gäller för små h (då h -> 0).
    % I detta fall avviker kurvan även då h blir för litet; det beror på
    % att referenslösningen inte är exakt och att vi hamnat för nära
    % referenslösningen.
end

function w = solve(h)
% Koden här är inklistrad från ENM_7_13.m.
w0 = [1 0 0 1]'; % initialvärden (kolonnvektor)
t0 = 0; tend = 10; % tidsintervall

% Högerled till ekvationen w' = f(w)
% Obs: f är en funktion från R^4 till R^4.
f = @(w) [w(2);
          -w(1)/(w(1)^2+w(3)^2)^(3/2);
          w(4);
          -w(3)/(w(1)^2+w(3)^2)^(3/2)];

N = (tend-t0) / h; % antal steg
t = t0:h:tend; % vektor med alla tidpunkter
w = zeros(length(w0),length(t)); % lösningsmatris
w(:,1) = w0; % sätt initialvärdet!
% Obs: varje kolonn i matrisen w motsvarar lösningen vid en viss tidpunkt.

% Euler framåt
for n=1:N
    w(:,n+1) = w(:,n) + h*f(w(:,n));
end
end
