% 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;