Metodi Numerici - Codici MATLAB svolti

Cap.11 Es.1

% Dato il problema u'(t)=1/(1+t^2)-2u^2 con u(0)=0 scrivere un codice matlab

% che implementa il metodo di Adams-Bashforth a 3 passi.

% Calcolare i primi passi in un caso con il metodo di Eulero e nell'altro

% con il metodo R-K esplicito del quarto ordine.

% Soluzione esatta u(t)=t/(1+t^2)

clear; clc; format long;

% Dati

t0=0; tf=3; u0=0; h=0.25;

t=t0;

fprintf(' t0 = %d, u0 = %d, h = %d \n ',t0,u0,h);

% Calcolo u1 con il metodi di Eulero

% Calcolo u2 con il metodo di R-K del quarto ordine (HH4)

% uab soluzione con Adams-Bashforth

% usol soluzione esatta

%calcolo del primo passo

f0=fval(t,u0);

n=1;

u1=u0+h*f0;

t=t+h;

usol=t/(1+t^2);

fprintf('\n n='); disp(n);

fprintf(' t='); disp(t);

fprintf(' u1 con Eulero='); disp(u1);

fprintf(' u(t1) sol. esatta='); disp(usol);

%calcolo del secondo passo

f1=fval(t,u1);

n=2;

u2=u1;

k1=fval(t,u2);

k2=fval(t+1/2*h,u2+1/2*h*k1);

k3=fval(t+1/2*h,u2+1/2*h*k2);

k4=fval(t+h,u2+h*k3);

u2=u2+(h/6)*(k1+2*k2+2*k3+k4);

t=t+h;

usol=t/(1+t^2);

fprintf('\n n='); disp(n);

fprintf(' t='); disp(t);

fprintf(' u1 con R-K del quarto ordine='); disp(u2);

fprintf(' u(t2) sol. esatta='); disp(usol);

%calcolo passi successivi

f2=fval(t, u2);

uab=u2;

while t<=tf-h

n=n+1

uab=uab+(h/12)*(23*f2-16*f1+5*f0);

%f2=f(tn,un);

%f1=f(tn-1,un-1);

%f0=f(tn-2,un-2);

f0=f1;

f1=f2;

f2=fval(t,uab);

usol=t/(1+t^2);

fprintf('\n n='); disp(n);

fprintf(' t='); disp(t);

fprintf(' u_Adams-Bashforth='); disp(uab);

fprintf(' sol. esatta='); disp(usol);

end

Cap.11 Es.2.2

% Scrivere i codici che implementano i metodi di Adams-Bashforth del secondo

% e del terzo ordine per il problema dell'esercizio 1.

clear; clc; format long;

% Dati

t0=0; tf=3; u0=0; h=0.25;

t(1)=t0;

fprintf(' t0 = %d, u0 = %d, h = %d \n ',t0,u0,h);

% Calcolo u(1) e u(2) con il metodo di Eulero

% uab2 soluzione con il metodo di Adams-Bashforth del secondo ordine

% uab3 soluzione con il metodo di Adams-Bashforth del terzo ordine

% usol soluzione esatta

%calcolo del primo passo

f0=fval(t(1),u0);

n=1;

u(1)=u0+h*f0;

t(2)=t(1)+h;

usol(1)=t(1)/(1+t(1)^2);

fprintf('\n n='); disp(n);

fprintf(' t='); disp(t(n));

fprintf(' u(1) con Eulero='); disp(u(1));

fprintf(' u(t1) sol. esatta='); disp(usol(1));

%calcolo del secondo passo

f1=fval(t(2),u(1));

X=f1; % Ci servirà poi per richiamare f1 in A-B secondo ordine

n=2;

u(2)=u(1)+h*f1;

t(3)=t(2)+h;

usol(2)=t(2)/(1+t(2)^2);

fprintf('\n n='); disp(n);

fprintf(' t='); disp(t(n));

fprintf(' u(2) con Eulero='); disp(u(2));

fprintf(' u(t2) sol. esatta='); disp(usol(2));

%calcolo passi successivi

f2=fval(t(3), u(2));

K=f2; % Ci servirà poi per richiamare f2 in A-B secondo ordine

n=3;

uab3(3)=u(2);

uab2(3)=u(2);

while t<=tf-h

n=n+1

uab3(n+1)=uab3(n)+(h/12)*(23*f2-16*f1+5*f0);

uab2(n+1)=uab3(n)+(h/2)*(3*K-X);

%f2=f(tn,un);

%f1=f(tn-1,un-1);

%f0=f(tn-2,un-2);

f0=f1;

f1=f2;

X=K;

t(n+1)=t(n)+h;

f2=fval(t(n+1),uab3(n));

K=fval(t(n+1),uab2(n));

usol(n+1)=t(n+1)/(1+t(n+1)^2);

fprintf('\n n='); disp(n);

fprintf(' t='); disp(t(n+1));

fprintf(' u_Adams-Bashforth3='); disp(uab3(n+1));

fprintf(' u_Adams-Bashforth2='); disp(uab2(n+1));

fprintf(' sol. esatta='); disp(usol(n+1));

end

hold on

plot(t,uab3,'ob');

plot(t,uab2,'og');

plot(t,usol,'or');

Cap.11 Es.2.3

% Scrivere i codici che implementano i metodi di Eulero implicito

% e dei Trapezi per il problema dell'esercizio 1.

clear; clc; format long;

% Dati

t0=0; tf=4; u0=0; h=0.1;

fprintf(' t0 = %d, u0 = %d, h = %d \n ',t0,u0,h);

% uimp soluzione con il metodo di Eulero implicito

% utrap soluzione con il metodo dei Trapezi

t(1)=t0; uimp(1)=u0; utrap(1)=u0

tol=1e-5; kmax=100; ymax=100; %tolleranza e numero massimo di iterazioni

fprintf( 'parametri del punto fisso \n');

fprintf( 'toll= %d, kmax=%d, ymax= %d, \n', tol, kmax, ymax);

n=1;

while t(n)<=tf-h

t(n+1)=t(n)+h;

fprintf('\n n='); disp(n);

fprintf(' t='); disp(t(n));

% Metodo Eulero implicito

v0=uimp(n); %inizializzazione dell'iterazione: u^(0)=u_n

% Ciclo k - Metodo del Punto Fisso

for k=1:kmax

v1=uimp(n)+h*(1/t(n+1)^2)-2*v0^2;

%criterio di arresto

if (abs(v1-v0)<tol*abs(v1));

break

end

v0=v1;

% variabile v0 assume valore u^(k)

% criterio di arresto non soddisfatto

end % fine ciclo Metodo del Punto Fisso

% Scrivere per controllo il numero iterazioni del Metodo del Punto Fisso

fprintf('iterazioni punto fisso per convergenza Eulero implicito='); disp (k);

uimp(n+1)=v1; % Calcolo di u_(n+1) valore di u^(k) che ha soddisfatto

% il criterio di arresto

% Metodo dei Trapezi

L0=utrap(n)

for y=1:ymax

L1=utrap(n)+(h/2)*((1/t(n)^2)-2*L0^2)+(1/(1+t(n+1)^2)-2*L0^2);

%criterio di arresto

if (abs(L1-L0)<tol*abs(L1));

break

end

L0=L1;

end

fprintf('iterazioni punto fisso per convergenza Trapezi='); disp (y);

utrap(n+1)=L1;

usol(n+1)=t(n+1)/(1+t(n+1)^2);

fprintf(' u_Eulero implicito='); disp(uimp(n+1));

fprintf(' u_Trapezi='); disp(utrap(n+1));

fprintf(' sol. esatta='); disp(usol(n+1));

end

hold on

plot(t,uimp,'ob');

plot(t,utrap,'og');

plot(t,usol,'or');

Cap.11 Es.2.1

% Scrivere i codici che implementano i metodi di Eulero, Eulero modificato,

% di R-K del quarto ordine per il problema dell'esercizio 1.

clear; clc; format long;

% Dati

t0=0; tf=3; u0=0; h=0.25;

t(1)=t0;

fprintf(' t0 = %d, u0 = %d, h = %d \n ',t0,u0,h);

% ueu soluzione con il metodo di Eulero

% urk soluzione con il metodo di R-K del quarto ordine (HH4)

% uem soluzione con il metodo di Eulero modificato

% usol soluzione esatta

ueu(1)=u0;

urk(1)=u0;

uem(1)=u0;

n=1;

while t(n)<=tf-h

ueu(n+1)=ueu(n)+h*fval(t(n),ueu(n));

k1=fval(t(n),urk(n));

k2=fval(t(n)+1/2*h,urk(n)+1/2*h*k1);