#!/usr/bin/octave
#testy otwartego Eulera dla r sklaranego y'=y y(0)=1 (rozw exp(t) 
# i y'=-y y(0)=1 (rozw exp(-t))
#oraz r 2 rzedu y''=-y z y(0)=1 y'(0)=0 (rozw cos(t))

clear all
t0=0;
n=10;T=1;h=(T-t0)/n;
x0=1; 
[x,t]=exEuler(@(x,t) x,x0,t0,h,n+1);
sol=@(x) exp(x);
plot(t,sol(t),";rozw;",t,x,";otw Euler;")
printf("y'=y [h=%g],ab=sol(%g)-x_n=%g, rel_err=%g\n",h,T,ar=sol(T)-x(n+1),ar/sol(T));
pause(3);

sol=@(x) exp(-x);
[x,t]=exEuler(@(x,t) -x,x0,t0,h,n+1);
plot(t,sol(t),";rozw;",t,x,";otw Euler;")
pause(3);

printf("y'=-y [h=%g],ab=sol(%g)-x_n=%g, rel_err=%g\n",h,T,ar=sol(T)-x(n+1),ar/sol(T));
#dluzszy czas

t0=0;n=100;T=10;
h=(T-t0)/n;
x0=1; 
[x,t]=exEuler(@(x,t) x,x0,t0,h,n+1);
sol=@(x) exp(x);
plot(t,sol(t),";rozw;",t,x,";otw Euler;")
printf("y'=y [h=%g],ab=sol(%g)-x_n=%g, rel_err=%g\n",h,T,ar=sol(T)-x(n+1),ar/sol(T));
pause(3);

sol=@(x) exp(-x);
[x,t]=exEuler(@(x,t) -x,x0,t0,h,n+1);
plot(t,sol(t),";rozw;",t,x,";otw Euler;")
pause(3);

printf("y'=-y [h=%g],ab=sol(%g)-x_n=%g, rel_err=%g\n",h,T,ar=sol(T)-x(n+1),ar/sol(T));



#testujemy otwartego Eulera dla y''=-y y(0)=1 y'(0)=0
function y=pw(x,t)
#pole wektorowe
 y=zeros(size(x));
 y(1)=x(2);
 y(2)=-x(1);
end


t0=0;
n=10;
T=1;h=(T-t0)/n;
x0=[1;0];
[x,t]=exEuler(@pw,x0,t0,h,n+1);
#wykres obl rozwiazan
plot(t,x(:,1),";y1;",t,x(:,2),";y2;");
pause(3);
plot(cos(t),-sin(t),";traj rozw;",x(:,1),x(:,2),";otw Eul t<1;");
pause(3);
sol=@(t) [cos(t);-sin(t)];
printf("y''=-y [h=%g],ab=||sol(%g)-x_n=%g||, rel_err=%g\n",\
        h,T,ar=norm(sol(T)-x(n+1,:)'),ar/norm(sol(T)));
T=10;n=100;h=(T-t0)/n;
 
[x,t]=exEuler(@pw,x0,t0,h,n+1);
#wykres obl rozwiazan
plot(t,x(:,1),";y1;",t,x(:,2),";y2;");
pause(3);
sol=@(t) [cos(t);-sin(t)];
plot(cos(t),-sin(t),";rozw;",x(:,1),x(:,2),";obl trajek t=[0,10];");
pause(3);

printf("y''=-y [h=%g],ab=||sol(%g)-x_n=%g||, rel_err=%g\n",\
        h,T,ar=norm(sol(T)-x(n+1,:)'),ar/norm(sol(T)));

#zmniejszamy h
clear x;
x0=1;
T=1;n=10*n;h=(T-t0)/n;
[x,t]=exEuler(@(x,t) x,x0,t0,h,n+1);
sol=@(x) exp(x);
plot(t,sol(t),";rozw;",t,x,";otw E h=1e-2;")
printf("y'=y [h=%g],ab=sol(%g)-x_n=%g, rel_err=%g\n",h,T,ar=sol(T)-x(n+1),ar/sol(T));
pause(3);

sol=@(x) exp(-x);
[x,t]=exEuler(@(x,t) -x,x0,t0,h,n+1);
sol=@(x) exp(-x);
plot(t,sol(t),";rozw;",t,x,";otw E h=1e-2;")
printf("y'=-y [h=%g],ab=sol(%g)-x_n=%g, rel_err=%g\n",h,T,ar=sol(T)-x(n+1),ar/sol(T));
pause(3);

#2 wymiary h=h/10
x0=[1;0];
[x,t]=exEuler(@pw,x0,t0,h,n+1);
#wykres obl rozwiazan
plot(t,x(:,1),";y1;",t,x(:,2),";y2;");
pause(3);
sol=@(t) [cos(t);-sin(t)];
plot(cos(t),-sin(t),";rozw;",x(:,1),x(:,2),";obl trajek t=[0,1];");
pause(3);
sol=@(t) [cos(t);-sin(t)];
printf("y''=-y [h=%g],ab=||sol(%g)-x_n=%g||, rel_err=%g\n",\
        h,T,ar=norm(sol(T)-x(n+1,:)'),ar/norm(sol(T)));

T=10;n=10*n;h=(T-t0)/n;
[x,t]=exEuler(@pw,x0,t0,h,n+1);
#wykres obl rozwiazan
plot(t,x(:,1),";y1;",t,x(:,2),";y2;");
pause(3);
sol=@(t) [cos(t);-sin(t)];
plot(cos(t),-sin(t),";rozw;",x(:,1),x(:,2),";otw E h=1e-2 t=[0,10];");
pause(3);
sol=@(t) [cos(t);-sin(t)];
printf("y''=-y [h=%g],ab=||sol(%g)-x_n=%g||, rel_err=%g\n",\
        h,T,ar=norm(sol(T)-x(n+1,:)'),ar/norm(sol(T)));
