%% class 15/12/2017
% pray and predators
% ex 17.0.10

clear all
close all
clc

a=4; b=2; c=3; d=3; 

f1=@(x,y,t)(a*x-b*x.*y);
f2=@(x,y,t)(d*x.*y-c*y);

x0=3; y0=3;

t0=0; tN=10; 

% in order to establish convergence w.r.t. N
N=(200:100:1000);
figure()
hold on
for i=1:length(N)

    [ x, y, t ] = RK2ndO_cpld1DE( t0, tN, N(i), x0, y0, f1, f2 );
    plot(t,x)
end
hold off



% we choose N=1000
N=1000;
[ x, y, t ] = RK2ndO_cpld1DE( t0, tN, N, x0, y0, f1, f2 );

figure()
subplot(1,2,1)
hold on
plot(t,x,'r')
plot(t,y,'b')
legend('x(t) pray','y(y) predators')
xlabel('t')
hold off

% critical points
cr1=[0 0];
cr2=[c/d a/b];

subplot(1,2,2)
hold on
plot(x,y,'r')
plot(cr1(1),cr1(2),'kx')
plot(cr2(1),cr2(2),'sg')
legend('orbits','unstable critical point','stable center')
quiver(x(1:50:N),y(1:50:N),f1(x(1:50:N),y(1:50:N),0),f2(x(1:50:N),y(1:50:N),0))
xlabel('x (prays)')
xlabel('y (predators)')
hold off

% let's consider more initial conditions so that to explore the entire
% phase diagram

x0=(0:0.6:2.5); y0=(0:0.6:2.5);
t0=0; tN=5; N=1000;

figure()
hold on
plot(cr1(1),cr1(2),'kx')
plot(cr2(1),cr2(2),'sg')
for i=1:length(x0)
    for j=1:length(y0)
        [ x, y, t ] = RK2ndO_cpld1DE( t0, tN, N, x0(i), y0(j), f1, f2 );
        plot(x,y,'r')
        quiver(x(1:50:N),y(1:50:N),f1(x(1:50:N),y(1:50:N),0),f2(x(1:50:N),y(1:50:N),0),0.5,'b')
    end
end
legend('unstable critical point','stable center','orbits')
xlabel('x (prays)')
xlabel('y (predators)')
axis([-0.5,4,-0.5,7])
hold off
        
 
%% revision of indefinite integrals with cum sum
clear all
close all
clc

N=50;
x=linspace(0,2*pi,N);
dx=x(2)-x(1)
g=sin(x);

il=(1:1:N-1); ir=(2:1:N);

f(1)=0;
f(ir)=cumsum( (g(il)+g(ir))/2 *dx );

figure()
hold on
plot(x,-cos(x)+1,'b')
plot(x,f,'rx')
legend('analytical','numerical')
xlabel('x')
hold off

%%
clear all
close all
clc

% if we want to create an external function that does an indefinite
% integral for any function g(s)

xi=0; xf=pi; N=100;
g=@(s)(exp(-s.^2).*cos(5*s))

[ f, x ] = ind_int_trpz( xi, xf, N, g )

figure()
hold on
plot(x,f,'r')
legend('indefinite integral of g(s)=e^{-s^2}cos(5s) ')
xlabel('x')
hold off

%% mean and variance
% exercise 18.1.1
clear all
close all
clc

x=[2 9 3 12 9];
y=[5 6 7 8 9];

x_mean=mean(x)
% equivalent to
sum(x)/length(x)

sigma_x=sqrt(var(x))
% equivalent to
sqrt( sum((x-x_mean).^2)/(length(x)-1) )

y_mean=mean(x)
sigma_y=sqrt(var(y))

%% ex 18.2.2
clear all
close all
clc


data=[9 2 8 7 1 3 4 8 3 2 3 10 9 3 8 9];

x=(1:1:10);
[p,x]=hist(data,x);

figure()
hold on
hist(data,x)
plot(x,p,'xr')
xlabel('x')
ylabel('occurrence')
hold off

figure()
hold on
bar(x,p)
plot(x,p,'xr')
xlabel('x')
ylabel('occurrence')
hold off

% if we normalise the histogram then we get the probability distribution
N=sum(p);

p=p/N;
figure()
hold on
bar(x,p)
plot(x,p,'xr')
xlabel('x')
ylabel('probability')
hold off

%% throwing a dice with 4 sides 10 times
clear all
close all
clc

N=10
data=randi(4,[1,N])

x=(1:1:4);

figure()
subplot(1,2,1)
[p x]=hist(data,x)
hold on
bar(x,p)
ylabel('occurrence (p(x) no normalised)')
hold off

prob=p/length(data)
subplot(1,2,2)
hold on
bar(x,prob)
ylabel('probability p(x)')
hold off
% check
sum(prob)
% gives 1

% if we throw it more times we get the uniform distribution
N=10000
data=randi(4,[1,N])

x=(1:1:4);
figure()
subplot(1,2,1)
[p x]=hist(data,x)
hold on
bar(x,p)
ylabel('occurrence (p(x) no normalised)')
hold off

prob=p/length(data)
subplot(1,2,2)
hold on
bar(x,prob)
ylabel('probability p(x)')
hold off
% check
sum(prob)
% gives 1
