% Jämförelse mellan felen i Runge-Kuttametoderna Eulers metod, Heuns metod 
% och fjärde ordningens Runge-Kutta 
% du/dt = -u^2 + t^2
% u(1) = 1


% Huvudprogram 
function [] = euler_heun_rungekutta()

disp('       |        Euler        |         Heun        |      Runge-Kutta')
disp('   h   |     Rh     E2h/Eh=2 |     Rh     E2h/Eh=4 |     Rh    E2h/Eh=16')

eR = []; % slutresultat för Euler
hR = []; % slutersultat för Heun
rkR = [];% slutresultat för Runge-Kutta
fig = 0;
for h = 0.5./(2.^(0:10))
	% Antal steg på intervallet [1 2]
	nsteg = (2 - 1)/h;

	fig = fig + 1;
	figure(fig)
	clf
	hold on
	title(['Steglängd ' num2str(h)], 'FontSize', 14)
	xlabel('t', 'FontSize', 14)
	ylabel('u', 'FontSize', 14)

	riktningsfalt(@f, [1 2], [1 2]);
	eu = euler(@f, [1 2], 1, nsteg);
	hu = heun(@f, [1 2], 1, nsteg);
	rku = rungekutta(@f, [1 2], 1, nsteg);
	
	% Approximera felet Eh med abs(R2h - Rh)
	eR = [eR, eu(end)];
	hR = [hR, hu(end)];
	rkR = [rkR, rku(end)];
	if (numel(eR) < 3)
		fprintf('%6.4f | %12.6e        | %12.6e        | %12.6e\n', ...
			h, eR(end), hR(end), rkR(end))
	else
		fprintf('%6.4f | %12.6e %6.3f | %12.6e %6.3f | %12.6e %6.3f\n', ...
			h, eR(end), abs(eR(end-2) - eR(end-1))/abs(eR(end-1) - eR(end)), ...
			hR(end), abs(hR(end-2) - hR(end-1))/abs(hR(end-1) - hR(end)), ...
			rkR(end), abs(rkR(end-2) - rkR(end-1))/abs(rkR(end-1) - rkR(end)))
	end
	pause
end
end

% Differentialekvationsfunktion
function dudt = f(t, u)

dudt = -u.^2 + t.^2;
end


% Rita riktningsfält och "sann" lösning
function true_u = riktningsfalt(fun, tgrans, ugrans)

% Rita riktningsfält
[tkoord, ukoord] = ...
	meshgrid(tgrans(1):0.05:tgrans(2), ugrans(1):0.05:ugrans(2));
pilt = ones(size(tkoord));
pilu = f(tkoord, ukoord);

quiver(tkoord, ukoord, pilt, pilu, 'k')
axis equal
axis tight
axis([tgrans ugrans])

% Rita "sann" lösning med ode45(f, [t_start t_end], y_start)
options = odeset('RelTol', 1e-13);
[true_t, true_u] = ode45(fun, tgrans, 1, options);
plot(true_t, true_u, 'k', 'LineWidth', 2)
end


% Eulers metod
function all_u = euler(fun, tgrans, ustart, nsteg)

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

% De t för vilka vi beräknar u
all_t = tgrans(1):h:tgrans(2);

% Ta steg h i t-led, beräkna u
t = tgrans(1);
u = ustart;
all_u(1) = u;
step = 1;
while t < tgrans(2)
	u = u + h*fun(t, u);
	t = t + h;
	
	step = step + 1;
	all_u(step) = u;
	
	plot(all_t([step-1 step]), all_u([step-1 step]), 'r', 'LineWidth', 2)
end
pause(0.1)
end


% Heuns metod
function all_u = heun(fun, tgrans, ustart, nsteg)

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

% De t för vilka vi beräknar u
all_t = tgrans(1):h:tgrans(2);

% Ta steg h i t-led, beräkna u
t = tgrans(1);
u = ustart;
all_u(1) = u;
step = 1;
while t < tgrans(2)
	k1 = h*fun(t, u);
	k2 = h*fun(t + h, u + k1);
	u = u + (k1 + k2)/2;
	t = t + h;
	
	step = step + 1;
	all_u(step) = u;
	
	plot(all_t([step-1 step]), all_u([step-1 step]), 'g', 'LineWidth', 2)
end
pause(0.1)
end


% Fjärde ordningens Runge-Kutta
function all_u = rungekutta(fun, tgrans, ustart, nsteg)

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

% De t för vilka vi beräknar u
all_t = tgrans(1):h:tgrans(2);

% Ta steg h i t-led, beräkna u
t = tgrans(1);
u = ustart;
all_u(1) = u;
step = 1;
while t < tgrans(2)
	k1 = h*fun(t, u);
	k2 = h*fun(t + h/2, u + k1/2);
	k3 = h*fun(t + h/2, u + k2/2);
	k4 = h*fun(t + h, u + k3);
	u = u + (k1 + 2*k2 + 2*k3 + k4)/6;
	t = t + h;
	
	step = step + 1;
	all_u(step) = u;
	
	plot(all_t([step-1 step]), all_u([step-1 step]), 'b', 'LineWidth', 2)
end
pause(0.1)
end

