%%%%%%%%%%%%%%%%%%%%%%
%
% Testowanie schematów 
%
%%%%%%%%%%%%%%%%%%%%%%


% definicja problemu 
%
clear;

f = @(t,x) x; % prawa strona równania, $f(t,x)$
t0 = 0; x0 = 1; T = 4; % dane zadania

% definicja rozwiązania dokładnego, $x(T)$ - dla sprawdzenia błędu
%
X = @(t) exp(t); 

% definicja parametrów schematu
%
N = 20; % liczba kroków 

%
% cz. I: test schematu
%
xe = euler(f, x0, t0, T, N);
xrk4 = rk4(f, x0, t0, T, N);

fprintf(stderr, 'Wartości w chwili T = %g\n\t Euler: %f | RK4: %f | dokł: %f\n', T, xe, xrk4, X(T)); % wartości w chwili T

fprintf(stderr, 'Błędy w chwili T = %g\n\t Euler: %e | RK4: %e | dokł: %e\n', T, [xe, xrk4,X(T)]-X(T)); % błędy w chwili T

%
% cz. II: wykres rozwiązania
%

t = linspace(t0, T, N+1); % punkty do wykresu
% wartości rozwiązania w chwilach $t_0,\ldots, T$ ...
xe = zeros(size(t)); xe(1) = x0; 	% ... dla Eulera
xrk4 = xe; 							% ... dla RK4

for i = 1:N
	% wykonujemy JEDEN krok schematu...
	xe(i+1) 	= euler(f, xe(i), t(i), t(i+1), 1 );
	xrk4(i+1) 	= rk4(f, xrk4(i), t(i), t(i+1), 1 );
end

fprintf(stderr, 'Wartości w chwili T = %g\n\t Euler: %f | RK4: %f | dokł: %f\n', T, xe(end), xrk4(end), X(T)); % wartości w chwili T

fprintf(stderr, 'Błędy w chwili T = %g\n\t Euler: %e | RK4: %e | dokł: %e\n', T, [xe(end), xrk4(end),X(T)]-X(T)); % błędy w chwili T

% rysujemy wykresy
plot(t,xe,'ro');
hold on;
plot(t,xrk4,'go');
plot(t, X(t), 'b-')
hold off;
legend('Euler', 'RK4', 'X(t)');

%
% cz. III: sprawdzenie rzędu
%
N = 10;
for i = 1:7 % uwaga: jeśli weźmiemy za dużo, może okazać się, że rząd 'nie wychodzi'
	xe = euler(f, x0, t0, T, N);
	xrk4 = rk4(f, x0, t0, T, N);
	err_e(i) = norm(xe-X(T));	
	err_rk4(i) = norm(xrk4-X(T));
	N = 2*N;
end

fprintf(stderr,'Oszacowanie rzędu schematu');
fprintf(stderr,'\n\tEulera:');
fprintf(stderr,' %f', log2(err_e(1:end-1)./err_e(2:end)));
fprintf(stderr,'\n\tRK4:   ');
fprintf(stderr,' %f', log2(err_rk4(1:end-1)./err_rk4(2:end)));
fprintf(stderr,'\n');
