% Questo programma risolve numericamente un problema di Cauchy, ossia un'equazione differenziale
% del tipo y'(t)=f(y(t)) con la condizione iniziale y(0)=y0. Per la risoluzione vengono utilizzati
% tre metodi:
% Metodo di Eulero esplicito. Approssimando la derivata con il rapporto incrementale si ottiene
% una semplice formula risolutiva di tipo iterativo u(i+1)=u(i)+f(u(i))*delta, dove delta indica
% la distanza tra due punti successivi tra quelli con cui e' stato suddiviso l'intervallo
% temporale in cui abbiamo deciso di lavorare.
% Metodo di Heun. Questo metodo, ottenuto come miglioramento del precedente, fornisce una formula
% del tipo u(i+1)=u(i)+1/2*[f(u(i))+f(u(i)+f(u(i))*delta)]*delta. La derivata nel punto corrente
% u(i) viene sostituita con la media aritmetica tra essa e quella nel punto corrispondente
% al passo di Eulero u(i)+f(u(i))*delta.
% Metodo di Eulero implicito. Stavolta la derivata viene approssimata tramite il rapporto
% incrementale da sinistra ottenendo u(i+1)=u(i)+f(u(i+1))*delta. In questo caso pero' otteniamo
% un'equazione del tipo u(i)+f(x)*delta-x=0, che viene qui risolta tramite il metodo di bisezione.




%                        Commenti a cura di Vittorio Maio




clear
close all
% Poiche' la funzione considerata in questo programma e' f(x)=c*x, cominciamo con il fissare y0 e c.
yo=1; 
c=1;
T=1; % [0,T] e' l'intervallo dove vogliamo approssimare la soluzione.

% Sappiamo che in questo caso la soluzione e' y(t)=y0*e^(c*t).


% Utilizzando un ciclo for (in j), possiamo far variare il numero n di intervallini di lunghezza
% delta tramite i quali abbiamo diviso l'intervallo [0,T]. In questo modo possiamo ottenere piu'
% risultati, di ciascuno dei quali faremo il plot insieme al plot dell'errore ottenuto ogni volta
% in funzione del passo delta. All'interno del primo ciclo for ne inseriamo un altro (in i) dove
% viene calcolata la soluzione approssimata tramite i tre metodi suddetti; dopo aver trovato la 
% soluzione numerica, la confrontiamo con la soluzione teorica a noi nota e calcoliamo l'errore
% massimo. Analizziamo poi la dipendenza di tale errore in funzione del passo delta. Poiche'
% ci aspettiamo una dipendenza di tipo legge di potenza, applichiamo il logaritmo ad entrambe
% le variabili per avere una relazione lineare, il cui coefficiente angolare ci fornisce l'ordine
% con cui va a zero l'errore.


for j=1:6
    n=10^j
    w=linspace(0,T,n+1);
    delta=1/(n);

    u=zeros(1,n+1);
    u(1)=yo;

    for i=1:n
%        u(i+1)=u(i)+c*u(i)*delta; %    Eulero
%       u(i+1)=u(i)+1/2*(c*u(i)+c*(u(i)+c*u(i)*delta))*delta; %    Heun
       u(i+1)=bisezione(c, delta , u(i)); %   Eulero implicito

    end


    z=yo*exp(c*w);
    figure
    plot(w,u,'.');
    hold on
    plot(w,z);

    numeri(j)=log10(delta);
    errori(j)=log10(max(abs(u-z)));
end


% Infine calcoliamo la retta dei minimi quadrati per ottenere l'ordine con cui va a zero l'errore
% tra la soluzione numerica e quella esatta in funzione di delta.


a=(errori-mean(errori))*(numeri-mean(numeri))'/((numeri-mean(numeri))*(numeri-mean(numeri))')
b=mean(errori)-a*mean(numeri)
z=a*numeri+b;
figure
plot(numeri,errori,'-')
hold on
plot(numeri,z,'+')


%--------------------------------------------------------------------------------------


%  mettere in un file chiamato bisezione.m



% Questa funzione MatLab riceve in input dal programma chiamante le variabili c, ui e delta ed in
% modo iterativo risolve nel seguente modo l'equazione g(x)=0, dove g(x):=ui+c*x*delta-x. Si parte
% da un intervallo iniziale [l,r] nei cui estremi la funzione g(x) ha segno discorde. La funzione
% g(x) viene calcolata nel punto medio (l+r)/2, dopodiche' si controlla se essa ha segno discorde
% rispetto ad g(l) o g(r). Nel primo caso si rimpiazza r con (l+r)/2, mentre nel secondo caso
% quest'ultimo punto sostituisce l. Iterando questo procedimento con un comando while si arriva
% alla radice dell'equazione con un errore inferiore ad una soglia prefissata eps. L'ouput e'
% la soluzione dell'equazione g(x)=0.


function [x] = bisezione(c,delta,ui)

eps=10^(-15);
l=ui;
r=2*ui;

while abs(r-l)>eps
    h=(l+r)/2;
    if (ui+c*h*delta-h)*(ui+c*l*delta-i)<0
        r=h;
    else l=h;
    end
end
x=h;