%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% numerical integration 
%
% example: int from 0 to 2 (x^2+1) dx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('exact: 14/3 = 4.6667')

N = 5;        % number of points (number of subintervals is N-1)
a = 0; b = 2; % endpoints of the given interval


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  solution using explicit loops
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h = (b-a)/(N-1);  % integration step

   % initialization of variables for the sums:
Sl = 0;    % Left Riemann    
Sr = 0;    % Right Riemann
Sm = 0;    % Middle-point rule
St = 0;    % Trapezoid rule

xl = a;             % left end of the interval
for k = 1:N-1       % loop for every subinterval
  xr = xl + h;      % right end of the subinterval
  xc = xl + h/2;    % the center of the subinterval
  Sl = Sl + h*(xl^2+1); 
  Sr = Sr + h*(xr^2+1);
  Sm = Sm + h*(xc^2+1);
  St = St + h*0.5*(xr^2+1 + xl^2+1);
  xl = xr;          % left end of the next subinterval
end

disp(['Left R. sum: ' , num2str(Sl)]);
disp(['Right R. sum: ' , num2str(Sr)]);
disp(['Midpoint rule: ' , num2str(Sm)]);
disp(['Trapezoid rule: ' , num2str(St) , ', ' , num2str(StM)]);


%%%%%%%%%%%%%%%%%%%%%%%%%
%  solution using vectors
%%%%%%%%%%%%%%%%%%%%%%%%%

x = linspace(a,b,N);  % equidistant points on x-axis in interval <a,b>

fx = x.^2+1;          % function values at given points

h = x(2)-x(1);        % step-size (or h = (b-a)/(N-1)  )

  % left Riemann sum
Sl = h*sum(fx(1:end-1));

  % right Riemann sum
Sr = h*sum(fx(2:end));

  % midpoint rule
xm = x + h/2;         % middle points
fm = xm.^2+1;
Sm = h*sum(fm(1:end-1));

  % trapezoid rule
ft = 0.5*(fx(1:end-1)+fx(2:end));
St = h*sum(ft);

  % trapezoid rule - Matlab built-in function
StM = h*trapz(fx);

disp(['Left R. sum: ' , num2str(Sl)]);
disp(['Right R. sum: ' , num2str(Sr)]);
disp(['Midpoint rule: ' , num2str(Sm)]);
disp(['Trapezoid rule: ' , num2str(St) , ', ' , num2str(StM)]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    logistic map
%    = a special case of the fixed point iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r = 3;     % parameter of the map
x = 0.5;   % initial size of the population (a number from 0 to 1)

  %%%  Problem 1: what will be the population size after N generations?

N = 10;    % number of generations
for k = 1:N
  x = r*x*(1-x);
end

disp(x)


  %%%  Problem 2: after how many generations a steady state occurs?

precision = 0.01;  % two states are considered equal with this precision
x_old = 0.5;       % initial size of the population

x_new = r*x_old*(1-x_old);
k = 1;             % the counter of the iterations

while abs(x_old-x_new) > precision
  x_old = x_new;
  x_new = r*x_old*(1-x_old);
  k = k+1;
end

disp(['Steady state is ', num2str(x_new), ', after ', num2str(k), ' iterations.' ])


  %%%  adding a safety-catch for infinite loop

r = 4.5;     % for this value no steady state is achieved

kmax = 10; % maximal number of iterations

precision = 0.01;  % two states are considered equal with this precision
x_old = 0.5;       % initial size of the population

x_new = r*x_old*(1-x_old);
k = 1;             % the counter of the iterations

while (abs(x_old-x_new) > precision) && k < kmax
  x_old = x_new;
  x_new = r*x_old*(1-x_old);
  k = k+1;
end

if k<kmax
  if isfinite(x_new)
    disp(['Steady state is ', num2str(x_new), ', after ', num2str(k), ' iterations.' ])
  else 
    disp('The process does not converge.');
  end
else 
  disp(['The process did not converge after ', num2str(k), ' steps.'])
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fixed point iterations for equation x = (x^3-7x+6)/7
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

precision = 0.0001;  % two values are considered equal with this precision
x_old = 0.5;        % initial guess

kmax = 3000; % maximal number of iterations

x_new = (x_old^3-7*x_old+6)/7;
k = 1;             % the counter of the iterations

while (abs(x_old-x_new) > precision) && k < kmax
  x_old = x_new;
  x_new = (x_old^3-7*x_old+6)/7;
  k = k+1;
end

if k<kmax
  if isfinite(x_new)
    disp(['A solution is ', num2str(x_new), ', after ', num2str(k), ' iterations.' ])
  else 
    disp('The process does not converge.');
  end
else 
  disp(['The process did not converge after ', num2str(k), ' iterations.'])
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%