% Jämförelse mellan Euler Framåt, Euler Bakåt, Trapetsmetoden ifråga om
% stabilitetsområde och noggrannhet

% Huvudprogram 
function [] = euler_trapets_stabilitet()

% GNM Exempel 6.5: Euler Bakåt, Trapetsmetoden stabil - Euler Framåt
% instabil. Trapetsmetoden också noggrannare i början!
% dy/dx = -8*x*y
% y(0) = 1
% Analytisk lösning y = e^(-4*x^2)
clear

% Definiera explicit funktion
f = @(x, y) -8*x*y;
% Definiera implicit funktion Euler Bakåt
g_eb = @(x, y, h) y/(1 + 8*h*(x + h)); 
% Definiera implicit funktion Trapetsmetoden
g_tr = @(x, y, h) y*(1 - 4*h*x)/(1 + 4*h*(x + h));

% Plotta analytisk lösning
true_x = 0:0.1:8;
true_y = exp(-4*true_x.^2);
 
figure(1)
clf
hold on
plot(true_x, true_y, 'k', 'LineWidth', 2)
xlabel('x')
ylabel('y')
axis([0 8 -1 1])
pause

% Antal steg
nsteg = 80;
% De x för vilka vi beräknar y
xgrans = [0 8];
% y(0)
ystart = 1;

% Euler Framåt
explicit_forstagrad(f, xgrans, ystart, nsteg, 'r');
pause
% Euler Bakåt
implicit_forstagrad(g_eb, xgrans, ystart, nsteg, 'g');
pause
% Trapetsmetoden
implicit_forstagrad(g_tr, xgrans, ystart, nsteg, 'b');
pause
axis([0 0.5 0.5 1])
end


% Euler Framåt
function all_y = explicit_forstagrad(fun, xgrans, ystart, nsteg, plotcol)

% Steglängd
h =	(xgrans(2) - xgrans(1))/nsteg;

% De x för vilka vi beräknar y
all_x = xgrans(1):h:xgrans(2);

% Ta steg h i x-led, beräkna y
x = xgrans(1);
y = ystart;
all_y(1) = y;
step = 1;
while x < xgrans(2) - 0.001
	y = y + h*fun(x, y);
	x = x + h;
	
	step = step + 1;
	all_y(step) = y;
	
	
	plot(all_x([step-1 step]), all_y([step-1 step]), plotcol, 'LineWidth', 2)
end
pause(0.1)
end


% Euler Bakåt och Trapetsmetoden
function all_y = implicit_forstagrad(implicitfun, xgrans, ystart, ...
	nsteg, plotcol)

% Steglängd
h =	(xgrans(2) - xgrans(1))/nsteg;

% De x för vilka vi beräknar y
all_x = xgrans(1):h:xgrans(2);

% Ta steg h i x-led, beräkna y
x = xgrans(1);
y = ystart;
all_y(1) = y;
step = 1;
while x < xgrans(2) - 0.001
	y = implicitfun(x, y, h);
	x = x + h;
	
	step = step + 1;
	all_y(step) = y;
	
	plot(all_x([step-1 step]), all_y([step-1 step]), plotcol, 'LineWidth', 2)
end
pause(0.1)
end
