
function esercizio5

% Oggetto: demo dei metodi di RK2
%
% Fissato "N", e posti "x0=0", "xfin=10", approssimare mediante il metodo RK2
% la soluzione "y" del problema di Cauchy:
%
% *  y'(x)=-y(x)+2*cos(x), x0 <= x <= xfin
% *  y(0) =1
%
% nei punti x(k)=x0+k*h, con h=(xfin-x0)/N, k=0,...,N.
%
% Si eseguano due esempi aventi i passi:
%
% 1. h=0.1
% 2. h=0.05
%
% e quali parametri del metodo RK si pongano
%
% gamma1=1/2, gamma2=1/2 alpha=1, beta=1.
%
% Noto che la soluzione e' "y(x)=sin(x)+cos(x)", calcolare nei due esempi
% gli errori assoluti compiuti nei punti x=2, x=4, x=6, x=8, x=10.
% Si tabulino gli errori richiesti in formato esponenziale, con una cifra
% prima della virgola e 2 dopo la virgola.
%
% Tratto da Atkinson, Han: "Elementary Numerical Analysis", p.414.

% Dati problema di Cauchy
x0=0; xfin=10;
y0=1; f=@(x,y) 2*cos(x)-y;

% Dati metodo
gamma1=1/2; gamma2=1/2; alpha=1; beta=1;

% Soluzione
y=@(x) sin(x)+cos(x);

% esempio 1.
N=100; h=(xfin-x0)/N;
errori=esempio(x0,y0,xfin,f,N,gamma1,gamma2,alpha,beta,y);

% esempio 2.
N=200; h=(xfin-x0)/N;
errori=esempio(x0,y0,xfin,f,N,gamma1,gamma2,alpha,beta,y);

% >> esercizio4
% 
%  	 x:  2 err: 1.93e-03
%  	 x:  4 err: 2.55e-03
%  	 x:  6 err: 5.81e-05
%  	 x:  8 err: 2.48e-03
%  	 x: 10 err: 2.13e-03
%  
% 
%  	 x:  2 err: 4.68e-04
%  	 x:  4 err: 6.25e-04
%  	 x:  6 err: 2.01e-05
%  	 x:  8 err: 6.04e-04
%  	 x: 10 err: 5.23e-04
%
% >>
 








function [x,u]=RK2(x0,y0,xfin,f,N,gamma1,gamma2,alpha,beta)

%--------------------------------------------------------------------------
% Oggetto:
% Implementazione dei metodi RK2
%--------------------------------------------------------------------------
% Input:
% x0: puto iniziale
% y0: valore della soluzione nel puto iniziale
% xfin: puto finale
% f: fuzione in due variabili, che descrive il problema di Cauchy:
%             * y'(x)=f(x,y(x)) x>=x0, x <= xfin
%             * y(x0)=y0
% N: numero di suddivisioni equisp. di [x0,xfin]
% gamma1,gamma2,alpha,beta: parametri RK2
%--------------------------------------------------------------------------
% Output:
% x: mesh di puti in cui e' approssimata la soluzione.
% u: approssimazione della soluzione del problema di Cauchy in "xn=x0+nh"
% con "n=0,...,N", dove "Nh=xfin-x0"
%--------------------------------------------------------------------------

x=linspace(x0,xfin,N+1);
u=[y0];

h=(xfin-x0)/N;

for n=2:N+1
    f0=feval(f,x(n-1),u(n-1));
    f1=feval(f,x(n-1)+alpha*h,u(n-1)+beta*h*f0);
    Fn=gamma1*f0+gamma2*f1;
    u(n)=u(n-1)+h*Fn;
end





function errori=esempio(x0,y0,xfin,f,N,gamma1,gamma2,alpha,beta,y)

[x,u]=RK2(x0,y0,xfin,f,N,gamma1,gamma2,alpha,beta);

v=2:2:10;
errori=[];

for k=v
    ii=find(x == k);
    if length(ii) > 0
        err=abs(y(k)-u(ii));
        errori=[errori; err];
        fprintf('\n \t x: %2.0f err: %1.2e',k,err);
    end
end

fprintf('\n \n');
