Funzioni MATLAB per la risoluzione di problemi differenziali di Cauchy

T=[];Y=[];

end

FUNZIONE METODO DI JACOBI

function [x,iter,residuo,rho] = Jacobi_ril(A,b,x0,omega,nmax,toll)

%

______________________________________________________________________________

%

% Function: Jacobi_ril

% _____________________

%

%

% Chiamata: [x,iter,residuo,rho] = Jacobi_ril(A,b,x0,omega,nmax,toll);

%

%

% La function 'Jacobi_ril' risolve un sistema lineare Ax = b utilizzando

% il metodo iterativo di Jacobi oppure il metodo di rilassamento in

% parallelo(JOR).

%

% Parametri di ingresso:

% ______________________

%

% A matrice non singolare del sistema, di dimensione nxn, con tutti

% gli elementi diagonali diversi da 0;

%

% b vettore colonna di lunghezza n contenente il termine noto del

% sistema;

%

% x0 vettore colonna di lunghezza n contenente l'approssimazione

% iniziale della soluzione del sistema con cui si innesca il

% metodo;

%

% omega parametro del metodo di rilassamento; quando si assume omega=1 si

% ottiene il metodo di Jacobi;

%

% nmax numero massimo di iterazioni imposte;

% toll tolleranza per il test d'arresto (sul residuo normalizzato);

%

%

% Parametri di uscita:

% ____________________

%

% x matrice (iter+1)xn contenente sulla (i+1)-esima riga

% l'approssimazione all'iterazione i-esima della soluzione

% del sistema lineare, i=0,1,...,iter;

%

% iter numero di iterazioni effettuate dal metodo per raggiungere la

% tolleranza desiderata sul residuo normalizzato;

%

% residuo vettore colonna di lunghezza iter+1, contenente, nella

% componente (i+1)-esima, la norma del residuo alla iterazione

% i-esima, i=0,1,...,iter;

%

% rho raggio spettrale della matrice B di iterazione del metodo.

%

%_____________________________________________________________________________

____

%

% Esempio di utilizzo della function applicata al sistema lineare Ax = b:

%

% A= 4 -1 0 0 0

% -1 4 -1 0 0

% 0 -1 4 -1 0

% 0 0 -1 4 -1

% 0 0 0 -1 4

%

% b= 3

% 2

% 2

% 2

% 3

%

% con soluzione esatta:

%

% alpha= 1

% 1

% 1

% 1

% 1

%

______________________________________________________________________________

_

% %

% Istruzioni Matlab

% clear all; clc

% n = 5;

% A = 4*eye(n)-diag(ones(n-1,1),1)-diag(ones(n-1,1),-1)

% b = [3; 2*ones(n-2,1) ; 3]

% x0 = zeros(n,1);

% omega = 1;

% nmax = 100; toll = 1e-12;

%

% % Chiamata function e risultati

% [x,iter,residuo,rho] = Jacobi_ril(A,b,x0,omega,nmax,toll);

%

% x =

% 0 0 0 0 0

% 0.7500 0.5000 0.5000 0.5000 0.7500

% 0.8750 0.8125 0.7500 0.8125 0.8750

% 0.9531 0.9063 0.9063 0.9063 0.9531

% 0.9766 0.9648 0.9531 0.9648 0.9766

% 0.9912 0.9824 0.9824 0.9824 0.9912

% 0.9956 0.9934 0.9912 0.9934 0.9956

% 0.9984 0.9967 0.9967 0.9967 0.9984

% 0.9992 0.9988 0.9984 0.9988 0.9992

% 0.9997 0.9994 0.9994 0.9994 0.9997

% 0.9998 0.9998 0.9997 0.9998 0.9998

% 0.9999 0.9999 0.9999 0.9999 0.9999

% 1.0000 1.0000 0.9999 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% ..............................................

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% % i puntini indicano iterazioni che non si riportano, il numero di

% % iterazioni effettuate è iter = 33

%_____________________________________________________________________________

___

[n,m]=size(A);

%___________________________________________________________________

D=diag(diag(A));

B=eye(n)-omega*inv(D)*A; %matrice di iterazione B del metodo

%___________________________________________________________________

rho=max(abs(eig(B)));

% Verifica della condizione necessaria e sufficiente di convergenza

if rho>=1

disp('*** Il metodo non converge ***')

iter=0; x=[]; residuo=[];

return

end

g=omega*inv(D)*b;

xv=x0;

normb=norm(b,inf);

if normb<=eps

normb=1;

end

res=norm(A*x0-b,inf)/normb;

residuo=res;

x=x0';i=1;

while res>toll & i<=nmax

xn=B*xv+g;

x=[x; xn'];

res=norm(A*xn-b,inf)/normb;

residuo=[residuo;res];

xv=xn;

i=i+1;

end

iter=i-1;

if i>nmax

disp('** il metodo non converge nel numero di iterazioni fissato **');

end

FUNZIONE METODO DI GAUSS SEIDEL

function [x,iter,residuo,rho] = Gauss_Seidel_ril(A,b,x0,omega,nmax,toll)

%_____________________________________________________________________________

_

%

% Function: Gauss_Seidel_ril

% __________________________

%

%

% Chiamata:[x,iter,residuo,rho] = Gauss_Seidel_ril(A,b,x0,omega,nmax,toll);

%

%

% La function 'Gauss_Seidel_ril' risolve un sistema lineare quadrato Ax = b

% utilizzando il metodo iterativo di Gauss-Seidel oppure il metodo di

% rilassamento in serie (SOR).

%

%

% Parametri di ingresso:

% ______________________

%

% A matrice del sistema non singolare, di dimensione nxn con tutti

% gli elementi diagonali diversi da zero;

%

% b vettore colonna di lunghezza n contenente il termine noto del

% sistema;

%

% x0 vettore colonna di lunghezza n contenente l'approssimazione

% iniziale della soluzione del sistema con cui si innesca il

% metodo;

%

% omega parametro di rilassamento (0 < omega < 2); quando si assume

% omega = 1 si ottiene il metodo di Gauss-Seidel;

%

% nmax numero massimo di iterazioni imposte;

%

% toll tolleranza per il test d'arresto (sul residuo normalizzato);

%

%

% Parametri di uscita:

% ____________________

%

% x matrice di dimensione(iter+1)xn contenente sulla (i+1)-esima

% riga l'approssimazione alla iterazione i-esima della soluzione

% del sistema lineare, i=0,1,...,iter;

%

% iter numero di iterazioni effettuate dal metodo per raggiungere la

% tolleranza desiderata sul residuo normalizzato;

%

% residuo vettore colonna di lunghezza iter+1, contenente, nella

% componente (i+1)-esima, la norma del residuo alla iterazione

% i-esima, i=0,1,...,iter;

%

% rho raggio spettrale della matrice B di iterazione del metodo.

%

%_____________________________________________________________________________

__

%

% Esempio di utilizzo della function applicata al sistema lineare Ax=b dove

%

% A= 4 -1 0 0 0

% -1 4 -1 0 0

% 0 -1 4 -1 0

% 0 0 -1 4 -1

% 0 0 0 -1 4

%

% b= 3

% 2

% 2

% 2

% 3

%

% con soluzione esatta:

%

% alpha= 1

% 1

% 1

% 1

% 1

%

______________________________________________________________________________

_

%

% % Istruzioni Matlab

% clear all; clc

% n = 5;

% A = 4*eye(n)-diag(ones(n-1,1),1)-diag(ones(n-1,1),-1)

% b = [3; 2*ones(n-2,1) ; 3]

% x0 = zeros(n,1);

% omega = 1;

% nmax = 100; toll = 1e-12;

%

% % Chiamata della function

%

% [x,iter,residuo,rho] = Gauss_Seidel_ril(A,b,x0,omega,nmax,toll);

%

% % Risultati

% x =

%

% 0 0 0 0 0

% 0.7500 0.6875 0.6719 0.6680 0.9170

% 0.9219 0.8984 0.8916 0.9521 0.9880

% 0.9746 0.9666 0.9797 0.9919 0.9980

% 0.9916 0.9928 0.9962 0.9985 0.9996

% 0.9982 0.9986 0.9993 0.9997 0.9999

% 0.9996 0.9997 0.9999 0.9999 1.0000

% 0.9999 0.9999 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

% 1.0000 1.0000 1.0000 1.0000 1.0000

%_____________________________________________________________________________

___

[n,m]=size(A);

%_____________________________________________________________________________

___

D=diag(diag(A));

OE=omega*tril(A,-1);

B=eye(n)-omega*inv(D+OE)*A; % matrice di iterazione del metodo

%_____________________________________________________________________________

___

rho=max(abs(eig(B)));

% Verifica della condizione necessaria e sufficiente di convergenza

if rho>=1

disp('*** Il metodo non converge ***')

iter=0; x=[]; residuo=[];

return

end

Au=D\triu(A,+1);

Al=D\tril(A,-1);

g=D\b;

xv=x0; xn=x0;

normb=norm(b,inf);

if normb<=eps

normb=1;

end

res=norm(A*x0-b,inf)/normb;

residuo=res;

x=x0'; i=1;

while res>toll & i<=nmax

for j=1:n

xn(j)=(1-omega)*xv(j)+omega*(g(j)-Al(j,:)*xn-Au(j,:)*xv);

end

res=norm(A*xn-b,inf)/normb;

x=[x; xn'];

residuo=[residuo;res];

xv=xn;

i=i+1;

end

iter=i-1;

if i>nmax

disp('** il metodo non converge nel numero di iterazioni fissato **')

end

FUNZIONE METODO DI RUNGE-KUTTA 4

function [T,Y] = Rungekutta4(t0,tmax,n,y0,f)

%

%___________________________________________________________________________

%

% Function Rungekutta4

% ____________________

%

% Chiamata: [T,Y] = Rungekutta4(t0,tmax,n,y0,f);

%

% Commento:

%

% La function Rungekutta4 utilizza il metodo di Runge-Kutta esplicito, a

% quattro stadi del quarto ordine , per determinare le approssimazioni

% della soluzione di un problema di Cauchy, nei nodi t(i),i = 1,2,...,n+1

% contenuti in un intervallo [t0,tmax]; t0 = t(1), tmax = t(n+1).

%

% Il problema (scalare o vettoriale) e' del tipo:

%

% y'(t) = f(t,y(t)) con t in (t0,tmax]

% y(t0) = yo

%

%

% dove:

% y(t) = (y1(t),y2(t),...,ym(t))

% f(t,y(t)) = (f1(t,y(t)),f2(t,y(t)),...,fm(t,y(t))

% y(t0) = (y1(t0),y2(t0),...,ym(t0)) = (yo1,yo2,...,yom) = yo

%

% con m = 1 nel caso scalare e m > 1 nel caso vettoriale.

%___________________________________________________________________________

%

%

% Parametri di input della function:

% ___________________________________

%

% t0,tmax estremi dell'intervallo di integrazione;

%

% n numero di sottointervalli [t(i),t(i+1)], i = 1,...,n in cui si

% desidera suddividere l'intervallo [t0,tmax];

%

% y0 vettore riga di lunghezza m contenente la condizione iniziale

% y(t0) = yo = [yo1,yo2,...,yom];

%

% f vettore colonna contenente le m stringhe f1, f2,..,fm

% contenenti le espressioni delle componenti di f scritte in

% funzione delle variabili t e [y(1),y(2),...,y(m)] corrispondenti

% a [y1(t),y2(t),...,ym(t)]; si puo' creare il vettore f utilizzando

% il comando'strvcat', nel seguente modo: f = strvcat(f1,f2,...,fm).

%

%

% Parametri di output della