%-----------------------------%
% Exercises unit 2 --- part 1 %
%-----------------------------%
%%
% Ex. 1.0.1
clear all
clc
M=[1 9 2 23; 2 8 13 2; -2 9 7 1; 45 2 82 -11]
M2=M;

for i=2:4
   M2(i,:)=M(i,:)-M(i-1,:); 
end
M2

%%
% alternatively a bit longer way of doing the same thing is
% performing a double loop for:
M3=M;
for i=2:4
    for j=1:4
    M3(i,j)=M(i,j)-M(i-1,j);
    end
end

M3

%%
%-----------------------------%
% Ex 1.0.2
clear all
clc
x=linspace(0,2,20);
for i=1:20
   f(i)=x(i)^2;
end
f

% compare with
f2=x.*x

%%
%-----------------------------%
% Ex 1.0.3
clear all
close all
clc
a=0; b=2*pi; N=100;
x=linspace(a,b,N);

for i=1:4
   figure(i)
   f=sin(i*x)        
   hold on
   plot(x,f,'-r')
   xlabel('x')
   ylabel('f')
   hold off
end

%%
%-----------------------------%
% Ex 1.0.3 - "Achilles and the Tortoise" plotted with loop "for"
close all
clear all
clc
n=(0:1:60);
% estimate of the distance at which A and T meet
dist_est=sum(20./(2.^n))
% estimate of the time after which A and T meet
time_est=sum(4./(2.^n))

% A's velocity [m/s] 
v_A=20/4
% T's velocity [m/s] 
v_T=10/4

% plot
a=0; b=9; N=200;
t=linspace(a,b,N);
fA=0+v_A*t;
fT=20+v_T*t;

% actual figure
hold on
plot(t,fA,'-b',...
    'LineWidth',3)
plot(t,fT,'-r',...
    'LineWidth',3)

one=ones(1,10);
for j=1:6
    if j==1
        tc(j,:)=linspace(0,4,10);
        plot(tc(1,:),20*one,'g','LineWidth',2)
        
        xc(j,:)=linspace(20,20+20/2,10);
        plot(4*one,xc(1,:),'-g','LineWidth',2)
    else
        n1=(0:1:j-2);
        n2=(0:1:j-1);
        n3=(0:1:j);
        
        tc(j,:)=linspace(sum(4./(2.^n1)),sum(4./(2.^n2)),10);
        plot(tc(j,:),sum(20./(2.^n2))*one,'g','LineWidth',2)
        
        xc(j,:)=linspace(sum(20./(2.^n2)),sum(20./(2.^n3)),10);
        plot(sum(4./(2.^n2))*one,xc(j,:),'-g','LineWidth',2)
    end
end

plot(time_est,dist_est,'sk','MarkerSize',14,'LineWidth',2)

xlabel('t [s]','fontsize', 16)
ylabel('x [m]','fontsize', 16)
title('Zeno s paradox of Achilles and the Tortoise','fontsize',16)
hold off

%%
%-----------------------------%
% Ex 1.0.4 -- evaluate the factorial of a given number
clear all
clc
% option 1
% initialization
f1=1
n=input('wait till i type the number n')
for i=1:n
    f1=f1*(n+1-i);
end
% the factorial of n is given by
f1

% equivalent option 2
f2=1;
for j=1:n
    f2=j*f2;
end
f2

% compare with built-in function
factorial(n)

%%
% this last option 3 was found by your class mate miguel:
% it is ingenious (because ans get's updated)
% but it is "dangerous":
clc
n=input('wait till i type the number n')
1;
for i=1:n
    ans*i; % carefull!!! you cannot define it as a variable
           % for example fact=ans*i
           % otherwise ans get's not updated to the last value
end
ans
%%
%-----------------------------%
% Ex 1.0.5
clear all
clc

A=[1 4; 2 5; 3 6]
B=[1 3 5; 2 4 6]
C1=A*B

% we want to redifine the same operation with the loop for:
for i=1:3
    for k=1:3
        % now we have to sum
        % initialise
        C2(i,k)=0;
        for j=1:2
            % sum repeated index j
            C2(i,k)=C2(i,k)+A(i,j)*B(j,k);
        end
    end
end
C2

D1=A.*B'
% or equivalently
BT=B';
for i=1:3
    for j=1:2
        D2(i,j)=A(i,j)*BT(i,j);
    end
end
D2
%%
%-----------------------------%
% Ex 1.0.6
clear all
close all
clc
% in this excercise we want to see how the numerical evaluation 
% of a geometric series tends asymptotically
% to its exact value

% Geometric series
% sum_{n=0}^{infinity} x^n = 1/(1-x)
x=1/2;
a_exact=1/(1-x);

% finite sums
N_st=1; N_end=20;
N=(N_st:1:N_end)

for i=1:length(N)
    n=(0:1:N(i));
    aN(i)=sum(1./2.^n);
end
aN
aN_exact=(1-x.^(N+1))/(1-x);

figure(1)
hold on
plot(N,aN,'-r','LineWidth',3)
plot(N,aN_exact,'og','MarkerSize',8,'LineWidth',2)
plot(N,a_exact*ones(length(N)),'--k')
legend('finite N numerical','finite N exact',...
        'asymptotic limit','Location', 'SouthEast')
axis([0 21 1.4 2.1])
xlabel('N','fontsize',16)
ylabel('aN','fontsize',16)
hold off

N_inv=1./N;
figure(2)
hold on
plot(N_inv,aN,'-r','LineWidth',3)
plot(N_inv,aN_exact,'og','MarkerSize',8,'LineWidth',2)
plot(N_inv,a_exact*ones(length(N)),'--k')
legend('finite N numerical','finite N exact',...
        'asymptotic limit','Location', 'SouthWest')
axis([0 1.02 1.48 2.02])
xlabel('N^{-1}','fontsize',16)
ylabel('aN','fontsize',16)
hold off

%%
%---------------------------%
% Example 1 on derivatives  %
%---------------------------%
% derivative of f(x)=x^2+5x at x0=2
clear all
close all
clc
x0=2;
deltax=0.1;
f2=(x0+deltax)^2+5*(x0+deltax);
f1=(x0)^2+5*(x0);

% approximated numerical value
fp_x0=(f2-f1)/deltax
% compare with exact value
fp_x0_exact=2*x0+5

%%
%------------------------%
% Exercise 2.0.7         %
%------------------------%
clear all
close all
clc

% derivative of f=x^2+5x at x0=2 for 
% different values of the increment

x0=2;

deltax=linspace(1.0e-6,1,100);
f2=(x0+deltax).^2+5*(x0+deltax);
f1=(x0)^2+5*(x0);
f3=(x0-deltax).^2+5*(x0-deltax);

fp_x0_left=(f2-f1)./deltax;
fp_x0_right=(f1-f3)./(deltax);
fp_x0_exact=(2*x0+5)*ones(1,length(deltax));

hold on
plot(deltax,fp_x0_left,'b-x')
plot(deltax,fp_x0_right,'r-o')
plot(deltax,fp_x0_exact,'k--')
%axis([1.0e-6 1 8.5 10.5])
xlabel('deltax','fontsize', 16)
ylabel('f`(x_0)','fontsize', 16)
hold off

%%
%---------------------------%
% Example 2 on derivatives  %
%---------------------------%
% f(x)=x^2+5x defined on interval [a,b]
clear all
close all
clc
a=0; b=10; N=30;
x=linspace(a ,b, N);
f=x.^2+5*x;
%figure(1)
%plot(x,f)

deltax=(b-a)/(N-1);
%deltax
%x(2)-x(1)

for i=1:N-1
    fp(i)=(f(i+1)-f(i))/(x(i+1)-x(i));
    xl(i)=x(i);
    xr(i)=x(i+1);
end
figure(2)
plot(xl,fp,'ro', xr,fp,'gx', x,2*x+5,'b')
legend(' approx RIGHT', ' approx LEFT', ' exact','Location','NorthWest')
xlabel('x'), ylabel('f')
%%
%equivalent to previous excercise
clear all
a=0; b=10; N=10;
x=linspace(a ,b, N);
f=x.^2+5*x;

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

%deltax(1)
%(b-a)/(N-1)

deltaf=f(ir)-f(il);

fp=deltaf./deltax;

plot(x(il),fp,'ro', x(ir),fp,'gx', x(il)+deltax/2,fp,'ks', x,2*x+5,'b')
legend(' approx 1', ' approx 2', ' approx 3', ' exact','Location','NorthWest')
xlabel('x'), ylabel('f')

%%
%------------------------%
% Exercise 2.0.8         %
%------------------------%
% f(x)=exp(x) defined on interval [a,b]
clear all
close all
clc
a=0; b=10; N=20;
x=linspace(a ,b, N);
f=exp(x);

deltax=(b-a)/(N-1);
%deltax
%x(2)-x(1)

for i=1:N-1
    fp(i)=(f(i+1)-f(i))/(x(i+1)-x(i));
    xl(i)=x(i);
    xr(i)=x(i+1);
end
figure(1)
plot(xl,fp,'ro', xr,fp,'gx', x,f,'b')
legend(' approx RIGHT', ' approx LEFT', ' exact','Location','NorthWest')
xlabel('x'), ylabel('f')
%%
%instead of introducing the index vectors ir and il, one can directly use
%the command diff
clear all
a=0; b=10; N=20;
x=linspace(a ,b, N);
f=exp(x);

deltax=diff(x);

deltaf=diff(f);

fp=deltaf./deltax;
il=(1:1:N-1); ir=(2:1:N);
figure(2)
plot(x(il),fp,'ro', x(ir),fp,'gx', x(il)+deltax/2,fp,'ks', x,exp(x),'b')
legend(' approx 1', ' approx 2', ' approx 3', ' exact','Location','NorthWest')
xlabel('x'), ylabel('f')

%%
%-----------------------------%
% Example 1
% higher approximation: derivative of f(x)=x^2+5x at x0=2
clear all
close all
clc
x0=2;
deltax=0.1;
f2=(x0+deltax)^2+5*(x0+deltax);
f1=(x0)^2+5*(x0);
f3=(x0-deltax)^2+5*(x0-deltax);

% approximated numerical value
fp_x0=(f2-f1)/deltax
% higher order: comment: why now is already exact!!
fp2_x0=(f2-f3)/(2*deltax)
% compare with exact value
fp_x0_exact=2*x0+5

%%
%------------------------%
% Exercise 3.0.10        %
%------------------------%
% another example (power higher than 2!)
% higher approximation: derivative of f(x)=x^3+5x^2 at x0=2
clear all
close all
clc
x0=2;
deltax=0.1;
f2=(x0+deltax)^3+5*(x0+deltax)^2;
f1=(x0)^3+5*(x0)^2;
f3=(x0-deltax)^3+5*(x0-deltax)^2;

% approximated numerical value
fp_x0_right=(f2-f1)/deltax
fp_x0_left=(f1-f3)/deltax
% higher order: comment: why now is already exact!!
fp2_x0=(f2-f3)/(2*deltax)
% compare with exact value
fp_x0_exact=3*x0^2+10*x0

% let's see how fast fp_x0 goes to the exact value fp_xo_exact
% when we are diminiscing the increment deltax

%deltax=linspace(1.0e-6,1,100);
deltax=(1.0e-6:0.01:1);
f2=(x0+deltax).^3+5*(x0+deltax).^2;
f1=(x0)^3+5*(x0)^2;
f3=(x0-deltax).^3+5*(x0-deltax).^2;

fp_x0_right=(f2-f1)./deltax;
fp_x0_left=(f1-f3)./deltax
fp2_x0=(f2-f3)./(2*deltax)
fp_x0_exact=(3*x0^2+10*x0)*ones(1,length(deltax));

hold on
plot(deltax,fp_x0_right,'b-x')
plot(deltax,fp_x0_left,'g-x')
plot(deltax,fp2_x0,'r-o')
plot(deltax,fp_x0_exact,'k--')
%axis([1.0e-6 1 31 44])
xlabel('deltax')
ylabel('f`(x_0)')
hold off

%%
%---------------------------------------%
% differentiate a non-periodic function %
%---------------------------------------%
%------------------%
% Exercise 3.0.11  %
%------------------%
% f(x)=x^2+5x^5 defined on interval [a,b]

clear all
close all
clc
a=0; b=10; N=20;
x=linspace(a ,b, N);
f=x.^2+5*x.^5;

deltax=(b-a)/(N-1);
%deltax
%x(2)-x(1)

for i=1:N-1
    fp(i)=(f(i+1)-f(i))/(x(i+1)-x(i));
    x1(i)=x(i);
    x2(i)=x(i+1);
    x3(i)=(x(i)+x(i+1))/2;
end
figure(1)
plot(x1,fp,'ro', x2,fp,'gx', x3,fp,'ks', x,2*x+25*x.^4,'b')
legend(' approx 1', ' approx 2', ' approx 3', ' exact','Location','NorthWest')
xlabel('x'), ylabel('f')
%%
%equivalent to previous excercise
clear all

a=0; b=10; N=20;
x=linspace(a ,b, N);
f=x.^2+5*x.^5;

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

%deltax(1)
%(b-a)/(N-1)

deltaf=f(ir)-f(il);

fp=deltaf./deltax;

figure(2)
plot(x(il),fp,'ro', x(ir),fp,'gx', x(il)+deltax/2,fp,'ks', x,2*x+25*x.^4,'b')
legend(' approx 1', ' approx 2', ' approx 3', ' exact','Location','NorthWest')
xlabel('x'), ylabel('f')

%%
% use of the command diff
clear all

a=0; b=10; N=20;
x=linspace(a ,b, N);
f=x.^2+5*x.^5;

dx=diff(x)
fp=diff(f)./dx

il=(1:1:N-1); ir=(2:1:N);
figure(3)
plot(x(il),fp,'ro', x(ir),fp,'gx', x(il)+dx(1)/2,fp,'ks', x,2*x+25*x.^4,'b')
legend(' approx 1', ' approx 2', ' approx 3', ' exact','Location','NorthWest')
xlabel('x'), ylabel('f')

%%
%----------------------------------%
% Exercise 4.0.13                  %
%----------------------------------%

clear all
close all
clc
% integrate the function f(x)=exp(x)
%between
a = 0; b = 1;
%using N points
N = 10;

x=linspace(a,b,N);
i1=(1:1:N-1);
i2=(2:1:N);

f=exp(x);

deltax=(b-a)/(N-1);
%check it is equal to
%x(2)-x(1)

I_l = sum( f(i1)*deltax );
I_r = sum( f(i2)*deltax );
I_c1 = sum( (( f(i1)+f(i2) )/2)*deltax );
I_c2 = sum( exp( (x(i1)+x(i2))/2 )*deltax );

I_l
I_r
I_c1
I_c2

% compare with the analytical solution
exp(b)-exp(a)

%compare the solutions in different formats
format long
I_l
I_r
I_c1
I_c2
exp(b)-exp(a)
%%
% use of trapz(x,y)
% help trapz

a = 0; b = 1; N = 10;
x=linspace(a,b,N);
f=exp(x);

trapz(x,f)
%%
%-------------------------------------%
% Same as above but with the loop for %
%-------------------------------------%
clear all
close all
clc
% integrate the function f(x)=exp(x)
%between
a = 0; b = 1;
%using N points
N = 100;

x=linspace(a,b,N);

f=exp(x);

deltax=(b-a)/(N-1);

I_l = 0;
I_r = 0;
I_c1 = 0;
I_c2 = 0;
for i=1:N-1
    I_l = I_l + f(i)*deltax;
    I_r = I_r + f(i+1)*deltax;
    I_c1 = I_c1 + (( f(i)+f(i+1) )/2)*deltax;
    I_c2 = I_c2 + exp( (x(i)+x(i+1))/2 )*deltax;
end
format long
I_l
I_r
I_c1
I_c2

% compare with the analytical solution
exp(b)-exp(a)

%%
%----------------------------------%
% Exercise 4.0.14                  %
%----------------------------------%
%Plot the above result as a function of N
clear all
a = 0; b = 1;

for N=10:100;

    x=linspace(a,b,N);
    i1=(1:1:N-1);
    i2=(2:1:N);

    f=exp(x);

    deltax=(b-a)/(N-1);

    I_l(N) = sum( f(i1)*deltax );
    I_r(N) = sum( f(i2)*deltax );
    I_c1(N) = sum( (( f(i1)+f(i2) )/2)*deltax );
    I_c2(N) = sum( exp( (x(i1)+x(i2))/2 )*deltax );
    exact(N)=exp(b)-exp(a);
end

n=(10:1:100);
hold on
plot(n,I_l(n),'rx');
plot(n,I_r(n),'g+');
plot(n,I_c1(n),'bs');
plot(n,I_c2(n),'co');
plot(n, exact(n),'k','Linewidth',1);
legend('approx l', ' approx r', ' approx c1', 'approx c2', 'exact','Location','NorthEast')
xlabel('x'), ylabel('f')
hold off

%%
%----------------------------------%
% Exercise 4.0.15                  %
%----------------------------------%
clear all
close all
clc

% plot the function f(x)=int_0^x ds exp(-s^2) cos(3s)
% for x in [0, 2*pi]

% value of the function at x=2*pi

a=0; b=2*pi; deltax=0.01;
x=(a:deltax:b);
N=length(x);
i1=(1:1:N-1);
i2=(2:1:N);

intgr=exp(-x.^2).*cos(3*x);

f_l  = sum( intgr(i1)*deltax )
f_r  = sum( intgr(i2)*deltax )
f_c1 = sum( (( intgr(i1)+intgr(i2) )/2)*deltax )
f_c2 = sum( ( exp(-((x(i1)+x(i2))/2).^2).*cos(3*(x(i1)+x(i2))/2) )*deltax )

%%
% now we repeat for different values of b!
close all
clear all
clc
af=0; bf=2*pi; Nf=100;
xf=linspace(af,bf,Nf);

f_l(1)=0; f_r(1)=0; f_c1(1)=0; f_c2(1)=0;

for i=2:Nf
    a=0; b=xf(i); deltax=0.01;
    x=(a:deltax:b);
    N=length(x);
    i1=(1:1:N-1);
    i2=(2:1:N);
    
    intgr=exp(-x.^2).*cos(3*x);
    
    f_l(i)  = sum( intgr(i1)*deltax );
    f_r(i)  = sum( intgr(i2)*deltax );
    f_c1(i) = sum( (( intgr(i1)+intgr(i2) )/2)*deltax );
    f_c2(i) = sum( ( exp(-((x(i1)+x(i2))/2).^2).*cos(3*(x(i1)+x(i2))/2) )*deltax );
end

figure(1)
hold on
plot(xf,f_l,'r--')
plot(xf,f_r,'g--')
plot(xf,f_c1,'k-x')
plot(xf,f_c2,'k-o')
Legend('approx left', 'approx right', 'approx c1', 'approx c2')
xlabel('x')
ylabel('f(x)')
hold off

%
% simple alternative with cumtrapz():
% thanks to one of your class mates!

intgr_xf=exp(-xf.^2).*cos(3*xf);
dx_xf=xf(2)-xf(1);
f_cumtrapz=cumtrapz(intgr_xf)*dx_xf;

figure(2)
hold on
plot(xf,f_l,'r--')
plot(xf,f_r,'g--')
plot(xf,f_c1,'k-x')
plot(xf,f_c2,'k-o')
plot(xf,f_cumtrapz,'x-r')
Legend('approx left', 'approx right', 'approx c1', 'approx c2','cumtrapz')
xlabel('x')
ylabel('f(x)')
hold off



%%
%----------------------------------%
% Exercise 4.0.16                  %
%----------------------------------%
clear all
close all
clc

% plot the function f(x)=int_0^x ds sin(s)/s
% for x in [0, 6*pi]

af=0; bf=6*pi; Nf=100;
xf=linspace(af,bf,Nf);

f_l(1)=0; f_r(1)=0; f_c1(1)=0; f_c2(1)=0;

for i=2:Nf
    a=0; b=xf(i); deltax=0.01;
    x=(a:deltax:b);
    N=length(x);
    i1=(1:1:N-1);
    i2=(2:1:N);
    
    for j=1:N
        if j==1
            intgr(1)=1.;
        else
            intgr(j)=sin(x(j))/x(j);
        end
    end
    f_l(i)  = sum( intgr(i1)*deltax );
    f_r(i)  = sum( intgr(i2)*deltax );
    f_c1(i) = sum( (( intgr(i1)+intgr(i2) )/2)*deltax );
    f_c2(i) = sum( ( sin(((x(i1)+x(i2))/2))./((x(i1)+x(i2))/2) )*deltax );
end

figure(1)
hold on
plot(xf,f_l,'r--')
plot(xf,f_r,'g--')
plot(xf,f_c1,'k-x')
plot(xf,f_c2,'k-o')
Legend('approx left', 'approx right', 'approx c1', 'approx c2')
xlabel('x')
ylabel('f(x)')
hold off