%% class 24/11/2017
%%
% Ex. 12.0.2 with Euler, modified Euler and 2nd order RK
clear all
close all
clc

% function defining the first derivative
% N.B. i define it as an anonymous function
% in general of t AND x
f=@(x,t)(-2*(t-5/2));

% initial and final conditions
t0=1;
x0=-9/4;
tN=4; N=5

% use Euler method for 1st order ODE
[ x_E, t ] = Euler_1ODE( t0, x0, tN, N, f );

x_exact=@(s)(-(s-5/2).^2);
% tangent at t2
xt2=x_exact(t(2));
x_tang=@(s)( xt2 +( -2*(t(2)-5/2) )*(s-t(2)) )

tplot=linspace(t0, tN, 100);

figure()
hold on
plot(tplot,x_exact(tplot),'r')
plot(t,x_E)
plot(tplot,x_tang(tplot),'r--')
legend('exact solution','numerical (1st order Euler)','tangent of exact sol at t=t(2)','Location','South')
plot(t(2),x_exact(t(2)),'rs')
axis([t0, tN, -2.5, 1.5])
hold off

% modified Euler
[ x_mE, t ] = modifiedEuler_1ODE( t0, x0, tN, N, f );

% 2nd order Runge-Kutta
[ x_2RK, t ] = RK2ndO_1ODE( t0, x0, tN, N, f );

figure()
hold on
plot(tplot,x_exact(tplot),'r')
plot(t,x_E)
plot(t,x_mE,'x-')
plot(t,x_2RK,'o--')
legend('exact solution','numerical (1st order Euler)',...
    'numerical (2nd order [modified] Euler))',...
    'numerical (2nd order Runge Kutta))','Location','South')
axis([t0, tN, -2.5, 1.5])
hold off
%%
% Ex. 12.0.3 with Euler, modified Euler and 2nd order RK
clear all
close all
clc

f=@(x,t)(-3*x);

% initial and final conditions
t0=0;
x0=1;
tN=1; N=7;

% use Euler method for 1st order ODE
[ x_E, t ] = Euler_1ODE( t0, x0, tN, N, f );
% modified Euler
[ x_mE, t ] = modifiedEuler_1ODE( t0, x0, tN, N, f );
% 2nd order Runge-Kutta
[ x_2RK, t ] = RK2ndO_1ODE( t0, x0, tN, N, f );

x_exact=@(s)( x0*exp(-3*s) );
tplot=linspace(t0, tN, 100);

figure()
hold on
plot(tplot,x_exact(tplot),'r','Linewidth',2)
plot(t,x_E)
plot(t,x_mE,'x-')
plot(t,x_2RK,'o--')
legend('exact solution','numerical (1st order Euler)',...
    'numerical (2nd order [modified] Euler))',...
    'numerical (2nd order Runge Kutta))','Location','NorthEast')
hold off

% error
err_E=max(abs(x_exact(t)-x_E))
err_mE=max(abs(x_exact(t)-x_mE))
err_2RK=max(abs(x_exact(t)-x_2RK))

% compare with
dt=t(2)-t(1);
dt
dt^2

%% same excercise as before but plotting the errors as a function of 1/N
clear all
close all
clc

f=@(x,t)(-3*x);

% initial and final conditions
t0=0;
x0=1;
tN=1; N=(15:1:100);

x_exact=@(s)( x0*exp(-3*s) );

for i=1:length(N)
    Ni=N(i);
    % use Euler method for 1st order ODE
    [ x_E, t ] = Euler_1ODE( t0, x0, tN, Ni, f );
    % modified Euler
    [ x_mE, t ] = modifiedEuler_1ODE( t0, x0, tN, Ni, f );
    % 2nd order Runge-Kutta
    [ x_2RK, t ] = RK2ndO_1ODE( t0, x0, tN, Ni, f );
    
    % error
    err_E(i)=max(abs(x_exact(t)-x_E));
    err_mE(i)=max(abs(x_exact(t)-x_mE));
    err_2RK(i)=max(abs(x_exact(t)-x_2RK));
    
    % to compare with
    dt(i)=t(2)-t(1);
    dt_sq(i)=dt(i)^2; 
end

% errors as a function of N
figure()
hold on
plot(N,err_E,'bx')
plot(N,err_mE,'rx')
plot(N,err_2RK,'ro')
plot(N,dt,'b--')
plot(N,dt_sq,'r--')
legend('error Euler', 'error modified Euler', 'error with RK', 'dt','dt^2')
xlabel('N')
hold off


% errors as a function of 1/N
figure()
hold on
plot(1./N,err_E,'bx')
plot(1./N,err_mE,'rx')
plot(1./N,err_2RK,'ro')
plot(1./N,dt,'b--')
plot(1./N,dt_sq,'r--')
legend('error Euler', 'error modified Euler', 'error with RK', 'dt','dt^2')
xlabel('1/N')
hold off

%% exercise 12.0.4 but with different initial conditions
% so that we can see that x0=0 is an UNSTABLE CRITICAL POINT
% while x0=5 is A STABLE CRITICAL POINT
close all
clear all
clc

f=@(x,t)(x.*(5-x));

% we have two equilibrium points
x01=0;
% is the unstable one
x02=5;
% is the stable one

figure()
y=linspace(-1,6,100);

% to plot the tangent
for i=1:100
    if f(y(i),0)<0
        a(i)=-1;
        b(i)=-(5-2*y(i));
    else
        a(i)=1;
        b(i)=(5-2*y(i));
    end
end

hold on
plot(y,f(y,0))
plot(x01,f(x01,0),'xb')
plot(x02,f(x02,0),'or')
legend('f(x)=dx/dt','unstable critical point','stable critical point',...
    'Location','South')
quiver(y(1:3:100),f(y(1:3:100),0),a(1:3:100),b(1:3:100),0.6)
plot(y,0*ones(1,100),'k--')
xlabel('x')
ylabel('f(x)=dx/dt')
hold off

x_ex=@(s,x1)( 5*x1*exp(5*s)./( 5-x1+x1*exp(5*s) ) );


% initial and final conditions
t0=0;
x0=[-1, 0, 1, 5, 6];
tN=3; N=100;

% for one initial condition we compare with exact solution
figure()
x0i=x0(3);
% 2nd order Runge-Kutta
[ x_2RK, t ] = RK2ndO_1ODE( t0, x0i, tN, N, f );
hold on
plot(t,x_ex(t,x0i),'b','Linewidth',2)
plot(t,x_2RK,'b--')
legend('exact x1=1','numerical RK x1=1')
xlabel('t')
ylabel('x(t)')
hold off

% then we look at the numerical one for a loop of initial conditions

figure()
hold on
for i=1:length(x0)
    
    x0i=x0(i);
    % 2nd order Runge-Kutta
    [ x_2RK, t ] = RK2ndO_1ODE( t0, x0i, tN, N, f );
    plot(t,x_2RK)
end
legend('RK x1=-1', 'RK x1=0', 'RK x1=1', 'RK x1=5', 'RK x1=6')
xlabel('t')
ylabel('x(t)')
axis([t0, tN, -3, 6])
hold off

