%----------------------%
% Example Euler method %
%----------------------%
%----------------%
% Exercise 2.0.1 %
%----------------%
clear all
close all
clc

t1=1; tM=4; M=5;
t=linspace(t1,tM,M);

dt=t(2)-t(1);
xp=-2*(t-5/2);

xt_E(1)=-9/4;

for i=1:M-1
   xt_E(i+1)=xt_E(i) + dt*xp(i);
end

% for plotting purposes only:
N=100;
tN=linspace(t1,tM,N);
xt_exact=-(tN-5/2).^2;

hold on
plot(tN,xt_exact,'r-','LineWidth',2)
plot(t,xt_E,'b-')
legend('exact solution', 'Euler approximation',...
    'Location','South')
xlabel('t','fontsize',20),ylabel('x(t)','fontsize',20)
hold off

print -f -depsc 'Euler_approx.eps'
%%
%----------------%
% Exercise 2.0.1 %
%----------------%
% Same as above, but using an anonymous functions

clear all
close all
clc

t1=1; tM=4; M=5;
% in general the first derivative f(x,t) can be a function
% of both x and t:
xp=@(x,t)( -2*(t-5/2) );

t=linspace(t1,tM,M);
dt=t(2)-t(1);

xt_E(1)=-9/4;

for i=1:M-1
   xt_E(i+1)=xt_E(i) + dt*xp(1,t(i));
end

xt_exact=@(t)( -(t-5/2).^2 );

N=100;
tN=linspace(t1,tM,N);

hold on
plot(tN,xt_exact(tN),'r-','LineWidth',2)
plot(t,xt_E,'b-')
legend('exact solution', 'Euler approximation',...
    'Location','South')
xlabel('t','fontsize',16),ylabel('x(t)','fontsize',16)
hold off

%%
%----------------%
% Exercise 2.0.1 %
%----------------%
% Same as above, but writing the Euler method part
% as an external function routine Euler_method.m
clear all
close all
clc

t1=1; tM=4; M=5;
x1=-9/4;
% in general the first derivative f(x,t) can be a function
% of both x and t and as such is also used in the 
% external function routine Euler_method.m
xp=@(x,t)( -2*(t-5/2) );
xt_exact=@(t)( -(t-5/2).^2 );

%call external function routine Euler_method.m
[xt_E, t] = Euler_method(t1, x1, tM, M, xp);

% used for plotting only
N=100;
tN=linspace(t1,tM,N);

hold on
plot(tN,xt_exact(tN),'r-','LineWidth',2)
plot(t,xt_E,'b-')
legend('exact solution', 'Euler approximation',...
    'Location','South')
xlabel('t','fontsize',16),ylabel('x(t)','fontsize',16)
hold off

%%
%----------------%
% Exercise 2.0.2 %
%----------------%
clear all
close all
clc

t1=0; tN=1; N=20;
x1=1;

t=linspace(t1,tN,N);
dt=t(2)-t(1);

xt_E(1)=x1;
for i=1:N-1
   xt_E(i+1)=xt_E(i) + dt*( -3*xt_E(i) );
end

hold on
plot(t,exp(-3*t),'r-','LineWidth',2)
plot(t,xt_E,'b-')
legend('exact solution', 'Euler approximation',...
    'Location','North')
xlabel('t','fontsize',16),ylabel('x(t)','fontsize',16)
hold off

%%
%----------------%
% Exercise 2.0.2 %
%----------------%
% Same as above, but now calling the external function
% routine Euler_method.m
clear all
close all
clc

t1=0; tN=1; N=20;
x1=1;
% rembember in this case too, even if now f(x,t) denpends
% on x only, in general can be a function of both
% and as such is defined and used in Euler_method.m
xp=@(x,t)( -3*x );

[xt_E, t] = Euler_method(t1, x1, tN, N, xp)

hold on
plot(t,exp(-3*t),'r-s','LineWidth',2)
plot(t,xt_E,'b-x')
legend('exact solution', 'Euler approximation',...
    'Location','North')
xlabel('t','fontsize',16),ylabel('x(t)','fontsize',16)
hold off

%%
%---------------------%
% Harmonic oscillator %
%---------------------%
%----------------%
% Exercise 3.0.3 %
%----------------%
clear all
close all
clc
a=0; b=4*pi; N=100;
t=linspace(a,b,N); dt=t(2)-t(1);
x0=1; v0=0;
x_E(1)=x0; v_E(1)=v0;

for i=1:N-1
   x_E(i+1)=x_E(i) + dt*v_E(i);
   v_E(i+1)=v_E(i) + dt*(-x_E(i));
end

delta=pi/2;
x_exact=sqrt(x0^2+v0^2)*sin(t+delta);

figure(1)
hold on
plot(t,x_exact,'r-','LineWidth',2)
plot(t,x_E)
legend('exact solution', 'Euler approximation',...
    'Location','NorthWest')
hold off

v_exact=sqrt(x0^2+v0^2)*cos(t+delta);

figure(2)
hold on
plot(x_exact,v_exact,'r-','LineWidth',2)
plot(x_E,v_E)
legend('exact solution', 'Euler approximation',...
    'Location','NorthWest')
hold off
%%
%----------------%
% Exercise 3.0.3 %
%----------------%
% Same as above but with an external function
% second_ord_Euler.m
clear all
close all
clc
a=0; b=4*pi; N=100;
t=linspace(a,b,N); dt=t(2)-t(1);
x0=1; v0=0;
dv=@(x,t)(-x);

[x_E, v_E, t] = second_ord_Euler(a, b, x0, v0, N, dv);

delta=pi/2;
x_exact=sqrt(x0^2+v0^2)*sin(t+delta);

figure(1)
hold on
plot(t,x_exact,'r-','LineWidth',2)
plot(t,x_E)
legend('exact solution', 'Euler approximation',...
    'Location','NorthWest')
hold off

v_exact=sqrt(x0^2+v0^2)*cos(t+delta);

figure(2)
hold on
plot(x_exact,v_exact,'r-','LineWidth',2)
plot(x_E,v_E)
legend('exact solution', 'Euler approximation',...
    'Location','NorthWest')
hold off

%%
%--------------------------------------------------%
% Excercise 3.1.1 for the problem 2.0.1            %
%--------------------------------------------------%
clear all
close all
clc
a=1; b=4; M=25;
f=@(x,t)( -2*(t-5/2) );
x0=-9/4;

xt_exact=@(t)( -(t-5/2).^2 );

[xt_E, t] = Euler_method(a, x0, b, M, f);

% accumulated error
abs(xt_E(M)-xt_exact(t(M)))
% compared with 
dt=t(2)-t(1);
(dt^2)*(M-1)


%%
% error
clear all
close all
clc
a=1; b=4;
xt_E(1)=-9/4;

for M=5:25;
    t=linspace(a,b,M);
    dt=t(2)-t(1);
    
    xt_exact=-(t-5/2).^2;
    xp=-2*(t-5/2);
    
    for i=1:M-1;
        xt_E(i+1)=xt_E(i) + dt*xp(i);
    end
    
    error1(M)=max( abs(xt_E-xt_exact) );
    deltat(M)=dt;
    error2(M)=((deltat(M))^2)*(M-1);
end

m=(5:1:M);
verror1=error1(m);
verror2=error2(m);
vdeltat=deltat(m);

figure(1)
hold on
plot(m,verror1,'r',m,verror2,'g',m,vdeltat,'b')
legend('max error','theoretical error','dt')
xlabel('M')
hold off

%convergence

%dM=diff(m);
%derror1=diff(verror1);
%derror2=diff(verror2);
%derror3=diff(vdeltat);

%conv1=derror1./dM;
%conv2=derror2./dM;
%conv3=derror3./dM;

%ml=(5:1:M-1); mr=(6:1:M);

%figure(2)
%hold on
%plot(ml+1/2,conv1,'r',ml+1/2,conv2,'g',ml+1/2,conv3,'b')
%legend('convergence 1','convergence 2','convergence 3')
%xlabel('M')
%hold off

