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