Metodi Numerici - Codici MATLAB svolti

Implementazione dei metodi di Eulero implicito e dei Trapezi

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));endhold onplot(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;% Datit0=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 Trapezit(1)=t0; uimp(1)=u0; utrap(1)=u0tol=1e-5; kmax=100; ymax=100; %tolleranza e numero massimo di iterazionifprintf( '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('<br>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)
    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));endhold onplot(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;% Datit0=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 esattaueu(1)=u0;urk(1)=u0;uem(1)=u0;n=1;while

t(n) <= tf - hueu(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);
k3 = fval(t(n)+1/2*h,urk(n)+1/2*h*k2);
k4 = fval(t(n)+h,urk(n)+h*k3);
urk(n+1) = urk(n) + (h/6)*(k1+2*k2+2*k3+k4);
uem(n+1) = uem(n) + h*(fval(t(n)+1/2*h,uem(n)+1/2*h*fval(t(n),uem(n)));
t(n+1) = t(n) + h;
usol(n+1) = t(n+1)/(1+t(n+1)^2);
fprintf('n='); disp(n);
fprintf('t='); disp(t(n+1));
fprintf('u_Eulero='); disp(ueu(n+1));
fprintf('u_R-K='); disp(urk(n+1));
fprintf('u_Eulero modificato='); disp(uem(n+1));
fprintf('sol. esatta='); disp(usol(n+1));
n = n+1;
end
hold on
plot(t,ueu,'ob');
plot(t,urk,'og');
plot(t,uem,'ok');
plot(t,usol,'or');

Cap.11 Es.3% Scrivere il codice del metodo predittore-correttore composto dal% metodo di Adams-Bashforth del terzo ordine (esplicito) come predittore e dal% metodo di Adams-Moulton del terzo ordine (implicito) come correttore% per il problema dell'esercizio 1.

clear;

clc; format long;
% Dati
t0=0; tf=10; u0=0; h=0.25;
t=t0;
fprintf('t0 = %d, u0 = %d, h = %d \n',t0,u0,h);

% Calcolo u1 con il metodo di Eulero
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 u2 con il metodo di Eulero
f1=fval(t,u1);
n=2;
u2=u1+h*f1;
t=t+h;
usol=t/(1+t^2);
fprintf('\n n=');
disp(n);
fprintf(' t=');
disp(t);
fprintf(' u2 con Eulero=');
disp(u2);
fprintf(' u(t2) sol. esatta=');
disp(usol);

% Applicazione del metodo
upi=u2;
while t<=tf-h
    up=upi+(h/12)*(23*f2-16*f1+5*f0);
    t=t+h;
    fp=fval(t,up); % Valore predetto
    uc=upi+(h/12)*(5*fp+8*f2-f1); % Correzione
    f3=fval(t,uc); % Valore corretto
    upi=uc;
    f0=f1;
    f1=f2;
    f2=f3;
    usol=t/(1+t^2);
    fprintf('\n n=');
    disp(n);

disp(n);
fprintf(' t=');
disp(t);
fprintf(' u_predittore=');
disp(up);
fprintf(' u_correttore=');
disp(uc);
fprintf(' sol. esatta=');
disp(usol);
n=n+1
end

Cap.11 Es.5% Consideriamo l'equazione u'=lambda*u soggetta alla condizione iniziale% u(0)=1. Si scrivano i codici Matlab che implementano il metodo di Eulero,% un metodo di R-K del secondo ordine e uno del quarto, il metodo di% Eulero implicito e il metodo dei Trapezi.

clear; clc; format long;

% Dati

t0=0; tf=3; u0=1; lambda=-1; h=0.02;
fprintf(' t0 = %d, u0 = %d, h = %d \n ',t0,u0,h);

% ueu soluzione con metodo di Eulero

% urk2 soluzione con metodo R-K del secondo ordine (Eulero migliorato)

% urk4 soluzione con metodo R-K del quarto ordine

% uimp soluzione con metodo di Eulero implicito

% utrap soluzione con metodo dei Trapezi

% usol soluzione esatta

t=t0; ueu=u0; urk2=u0; urk4=u0; uimp=u0, utrap=u0;
n=0;
while t<=tf-h
n=n+1
% Eulero esplicito
eu=ueu+h*(lambda*ueu)
% R-K

2k1=lambda*urk2k2=lambda+(urk2+h*k1)urk2=urk2+(h/2)+(k1+k2)% R-K 4k1=lambda*urk4k2=lambda+(urk4+(h/2)*k1)k3=lambda+(urk4+(h/2)*k2)k4=lambda+(urk4+h*k3)urk4=urk4+(h/6)*(k1+2*k2+2*k3+k4);% Eulero implicitouimp=1/(1-h*lambda)*uimp% Trapeziutrap=(1+(lambda*h)/2)/(1-(lambda*h)/2)*utrap;t=t+h;usol=exp(lambda*t);fprintf('n='); disp(n);fprintf(' t='); disp(t);fprintf(' u_Eulero='); disp(ueu);fprintf(' u_R-K 2='); disp(urk2);fprintf(' u_R-K 4='); disp(urk4);fprintf(' u_Eulero implicito='); disp(uimp);fprintf(' u_Trapezi='); disp(utrap);fprintf(' sol. esatta='); disp(usol);endCap.13 Es.4% Considero l'equazione u'=-50*(u-cos(t)) soggetta alla condizione iniziale u(0)=0% Scrivere i codici Matlab per il metodo di Eulero e Eulero implicitoclear; clc; format long;u0=0; t0=0; tf=10; h=0.001;% ueu soluzione calcolata con il metodo di Eulero% uimp soluzione calcoata con il metodo di Eulero implicitot=t0;

ueu=u0; uimp=u0;
n=0;
while t<=tf-h
    n=n+1;
    t=t+h;
    ueu=ueu+h*(-50*(ueu-cos(t)));
    uimp=(1/(1+50*h))*(uimp+50*h*cos(t));
    fprintf('<br> n=');
    disp(n);
    fprintf(' u_Eulero= ');
    disp(ueu);
    fprintf(' u_Eulero implicito= ');
    disp(uimp);
end

Cap.13 Es.5% Scrivere codici Matlab che implementano metodi alle differenze finite% a scelta per la risoluzione del problema.

clear; clc; format long;
t0=0; u0=[1; 0]; tf=0.5; h=0.01;
A=[-16 12; 12 -9];
[q, d]=eig(A);
fprintf(' matrice con diagonale autovalori= ');
disp(d);
[q, d]=eig(A);
fprintf(' matrice con diagonale autovalori= ');
disp(d);
t=t0;
b=[16*cos(t)-13*sin(t); -11*cos(t)+9*sin(t)];
ueu=u0; uimp=u0;
n=0;
while t<=tf-h
    n=n+1;
    t=t+h;
    ueu=ueu+h*(A*u+b);
    uimp=(eye(2)-h*A)/uimp+h*b;
    d=uimp-ueu;
    fprintf('<br> d=');
    disp(d);
end

` per mantenere la formattazione del testo come nel codice originale e il tag `
` per inserire un paragrafo di testo. Inoltre, utilizza l'entità HTML `<` per rappresentare il carattere `<` e l'entità HTML `>` per rappresentare il carattere `>`.
n='); disp(n);fprintf(' u_Eulero= '); disp(ueu);fprintf(' u_Eulero implicito= '); disp(uimp);fprintf(' differenza= '); disp(d);endCap.13 Es.6% Scrivere codici Matlab che implementano metodi alle differenze finite% a scelta per la risoluzione del problema.clear; clc; format long;t0=0; u0=[1; 0; -1]; tf=0.5; h=0.001;A=[-21 19 -20; 19 -21 20 ; 40 -40 -40];% Calcolo degli autovalori della matrice del sistema[q, d]=eig(A);fprintf(' matrice con diagonale autovalori= '); disp(d);% Siccome l'autovalore è lambda=-16, la scelta del passo h% deve rispettare la condizione h*|lambda|<2t=t0;% ueu soluzione calcolata con Eulero esplicito (NON A-stabile)% uimp soluzione calcolata con Eulero implicito (A-stabile)% usol soluzione esatta;ueu=u0; uimp=u0;n=0;while
t<=tf-hn=n+1;t=t+h;ueu=ueu+h*A*u;uimp=(eye(2)-h*A)/uimp;usol=[+(1/2)*exp(-2*t)+(1/2)*exp(-40*t)*(cos(40*t)+sin(40*t));+(1/2)*exp(-2*t)-(1/2)*exp(-40*t)*(cos(40*t)+sin(40*t));-exp(-40*t)*(cos(40*t)+sin(40*t))]fprintf('\n n='); disp(n);fprintf(' u_Eulero= '); disp(ueu);fprintf(' u_Eulero implicito= '); disp(uimp);fprintf(' sol. esatta= '); disp(usol);endCap.13 Es.10% Implementare in Matlab il metodo di Rosenbrock% con la formula di Rosenbrockper il problema dell'es.5clear; clc; format long;t0=0; u0=[1; 0]; tf=0.5; h=0.01;A=[-16 12; 12 -9];% Calcolo degli aitovalori della matrice del sistema[q, d]=eig(A);fprintf(' matrice con diagonale autovalori= '); disp(d);% Siccome l'autovalore è lambda=-16, la scelta del passo h% deve rispettare la condizione h*|lambda|<2% Coefficienti di Rosenbrockalpha1=1+sqrt(6)/6; alpha2=1-sqrt(6)/6;b21=(-6-sqrt(6)+sqrt(58+20*sqrt(6)))/(6+2*sqrt(2));
beta21=b21;c2=(sqrt(6)/6+1/2)/(2*sqrt(6)/6-b21); c1=1-c2;t=t0;b=[16*cos(t)-13*sin(t); -11*cos(t)+9*sin(t)];% ueu soluzione calcolata con Eulero esplicito (NON A-stabile)% uimp soluzione calcolata con Eulero implicito (A-stabile)% urb soluzione calcolata con il metodo di Rosenbrockueu=u0; uimp=u0; urb=u0;n=0;while t<=tf-hn=n+1;%calcolo k1f1=A*urb+b;