Codici Matlab

FORMULA DEI TRAPEZI SEMPLICE E COMPOSTA – MAIN

% calcolo l'integrale di cos(x) tra 0 e pi/2

a = 0;

b = pi/2;

n = 10;

fun = @(x) cos(x);

I_T = trapezi_semplice(fun,a,b);

I_V = sin(b)-sin(a);

I_Tn = trapezi_composta(fun,a,b,n);

I_Tn2 = trapezi_composta2(fun,a,b,n);

fprintf('integrale con trapezi: %f\n', I_T);

fprintf('integrale vero %f\n', I_V);

fprintf('integrale trapezi composto: %f\n', I_Tn);

fprintf('integrale trapezi composto2: %f\n', I_Tn2)

TRAPEZI SEMPLICE – FUNCTION

function integrale = trapezi_semplice(fun,a,b)

integrale = ((b-a)/2)*(fun(a)+fun(b));

end

TRAPEZI COMPOSTA – FUNCTION (METODO 1)

function integrale = trapezi_composta(fun,a,b,n)

% definisco il sottointervallo

h = (b-a)/n;

% inizializzo il ciclo

integrale = 0;

xi = a;

for i = 1:n

xip1 = xi+h;

Is = trapezi_semplice(fun, xi, xip1);

integrale = integrale + Is;

xi = xip1;

end

end

TRAPEZI COMPOSTA – FUNCTION (METODO 2)

function integrale = trapezi_composta2(fun,a,b,n)

h = (b-a)/n;

x = a:h:b;

f = fun(x);

integrale = h*( (f(1) + f(n+1) )/2 + sum(f(2:n)));

end

CAVALIERI SIMPSON COMPOSTA + ESTRAPOLAZIONE DI RICHARDSON – MAIN

clear

clc

a = 1;

b = 5;

fun = @(x) sin(x) + x.^2;

FUN = @(x) -cos(x) + (x.^3)/3;

% Calcolo integrale "vero"

IV = FUN(b) - FUN(a);

% Calcolo integrali CS

CS_10 = CavSim(10,a,b,fun);

CS_20 = CavSim(20,a,b,fun);

CS_40 = CavSim(40,a,b,fun);

% Calcolo con estrapolazione

R_10 = EstRic(CS_10,CS_20);

R_20 = EstRic(CS_20,CS_40);

% Calcolo errori

E_CS_10 = abs(IV-CS_10);

E_CS_20 = abs(IV-CS_20);

E_CS_40 = abs(IV-CS_40);

E_R_10 = abs(IV-R_10);

E_R_20 = abs(IV-R_20);

% Stampa risultati

fprintf('CS_10: %15.6e | Errore: %15.6e\n',CS_10,E_CS_10);

fprintf('CS_20: %15.6e | Errore: %15.6e\n',CS_20,E_CS_20);

fprintf('CS_40: %15.6e | Errore: %15.6e\n',CS_40,E_CS_40);

fprintf('R_10: %15.6e | Errore: %15.6e\n',R_10,E_R_10);

fprintf('R_20: %15.6e | Errore: %15.6e\n',R_20,E_R_20);

return

CAVALIERI SIMPSON - FUNCTION

function I = CavSim(n,a,b,fun)

h = (b-a) / n;

x = linspace(a,b,2*n+1);

fx = fun(x);

I = (h/6)*(fx(1)+fx(end)+4*sum(fx(2:2:end-1))+2*sum(fx(3:2:end-2)));

end

ESTRAPOLAZIONE DI RICHARDSON - FUNCTION

function I = EstRic(Q1,Q2)

I = Q2 + (Q2-Q1)/15;

end

CAVALIERI SIMPSON - MAIN

clear

close all

clc

% Definiamo estremi integrazione

a = 1;

b = 2;

% Definiamo funzione integranda

fun = @(x) x.*log(x);

% Definisco primitiva

% int[x*log(x)] = x^2/2*log(x) - int[x^2/2 * 1/x] = (x^2/2)*log(x) - x^2/4

FUN = @(x) 0.5*(x.^2).*log(x) - x.^2/4;

% Calcolo integrale vero

Ivero = FUN(b) - FUN(a);

% Faccio variare i il numero di sottointervalli

IQ = zeros(8,1);

EQ = zeros(8,1);

EQ_old = 1;

fprintf('nsub Ivero IQ EQ EQ_old/EQ\n');

for pow = 0:7

% Calcolo il numero di sottointervalli

nsub = 2^pow;

IQ(pow+1) = cavsimcomp(nsub,a,b,fun);

EQ(pow+1) = abs(IQ(pow+1)-Ivero);

fprintf('%d %f %f %e %f\n',...

nsub,Ivero,IQ(pow+1),EQ(pow+1),EQ_old/EQ(pow+1));

EQ_old = EQ(pow+1);

end

semilogy(EQ);

CAVALIERI SIMPSON - FUNCTION

function I = cavsimcomp(nsub,a,b,fun)

h = (b-a)/nsub;

x = a:h/2:b;

y = fun(x);

sum_medi = sum(y(2:2:end-1));

sum_sepa = sum(y(3:2:end-2));

I = (h/6)*(y(1)+y(end)+4*sum_medi+2*sum_sepa);

end

JACOBI – MAIN

clear

clc

% leggi soluzione e termine noto

A = load('matrice.dat');

b = load('tnoto.dat');

n = size(A,1);

x0= zeros(n,1);

toll = 1.e-8;

itmax = 150;

fout = fopen('risul_jacobi.dat', 'w');

% chiamata a jacobi

[xnew, iter, scarto, asint, R] = jacobi_comp(A, b, x0, toll,

itmax);

[xnew, iter, scarto, asint, R] = jacobi_mat(A, b, x0, toll,

itmax);

semilogy(2:iter, v_scarti(2:iter));

% calcolo autovalori

D_inv = inv(diag(diag(A)));

EJ = eye(n)-D_inv*A;

r_spett = max(abs(eig(EJ)));

JACOBI_MAT – FUNCTION

function [xnew, iter, scarto, asint, R] = jacobi_mat(A, b, x0,

toll, itmax)

n = size(A,1);

x_old = x0;

v_scarti = zeros(itmax,1);

% calcolo EJ

EJ = eye(n)-inv(diag(diag(A)))*A;

qj = inv(diag(diag(A)))*b;

iter = 0;

scarto = 2.0*toll;

scarto_old = 1.0;

while scarto > toll && iter < itmax

iter = iter + 1;

xnew = EJ*xold+qJ;

scarto = norm(xnew-xold);

asint = scarto/scarto_old;

% memorizzo lo scarto nel vettore degli scarti

v_scarti(iter) = scarto;

% passo all'iterazione sucessiva

xold = xnew;

scarto_old = scarto;

end

R = -log10(asint);

end

JACOBI_COMP – FUNCTION

function [xnew, iter, scarto, asint, R] = jacobi_comp(A, b, x0,

toll, itmax)

n = size(A,1);

x_old = x0;

v_scarti = zeros(itmax,1);

iter = 0;

scarto = 2.0*toll;

scarto_old = 1.0;

while scarto > toll && iter < itmax

iter = iter + 1;

% ciclo sulle componenti di xnew

for i = 1:n

som_L = 0.0;

for j = 1:i-1

som_L = som_L + A(i,j)*x_old(j);

end

som_U = 0.0;

for j = i+1:n

som_U = som_U + A(i,j)*xold(j);

end

xnew(i) = ((b(i)-som_L-som_U)/A(i,i));

end

scarto = norm(xnew-xold);

asint = scarto/scarto_old;

% memorizzo lo scarto nel vettore degli scarti

v_scarti(iter) = scarto;

% passo all'iterazione sucessiva

xold = xnew;

scarto_old = scarto;

end

R = -log10(asint);

end

SEIDEL – MAIN

clear

clc

% leggi soluzione e termine noto

A = load('matrice.dat');

b = load('tnoto.dat');

n = size(A,1);

x0= zeros(n,1);

toll = 1.e-8;

itmax = 150;

fout = fopen('risul_jacobi.dat', 'w');

% chiamata a seidel

[xnew, iter, scarto, asint, R] = seidl_com(A, b, x0, toll, itmax);

[xnew, iter, scarto, asint, R] = seidel_mat(A, b, x0, toll,

itmax);

semilogy(2:iter, v_scarti(2:iter));

% calcolo autovalori

L = tril(A);

ES = eye(n)-inv(L)*A;

r_spett = max(abs(eig(ES)));

SEIDEL_MAT – FUNCTION

function [xnew, iter, scarto, asint, R] = seidel_mat(A, b, x0,

toll, itmax)

n = size(A,1);

x_old = x0;

v_scarti = zeros(itmax,1);

% calcolo ES

L = tril(A);

ES = eye(n)-inv(L)*A;

qS = inv(L)*b

iter = 0;

scarto = 2.0*toll;

scarto_old = 1.0;

while scarto > toll && iter < itmax

iter = iter + 1;

xnew = ES*xold+qS;

scarto = norm(xnew-xold);

asint = scarto/scarto_old;

% memorizzo lo scarto nel vettore degli scarti

v_scarti(iter) = scarto;

% passo all'iterazione sucessiva

xold = xnew;

scarto_old = scarto;

end

R = -log10(asint);

end

SEIDEL_COMP – FUNCTION

function [xnew, iter, scarto, asint, R] = seidel_comp(A, b, x0,

toll, itmax)

n = size(A,1);

x_old = x0;

v_scarti = zeros(itmax,1);

iter = 0;

scarto = 2.0*toll;

scarto_old = 1.0;

while scarto > toll && iter < itmax

iter = iter + 1;

% ciclo sulle componenti di xnew

for i = 1:n

som_L = 0.0;

for j = 1:i-1

som_L = som_L + A(i,j)*x_old(j);

end

som_U = 0.0;

for j = i+1:n

som_U = som_U + A(i,j)*xold(j);

end

xnew(i) = (b(i)-som_L-som_U)/A(i,i);

end

scarto = norm(xnew-xold);

asint = scarto/scarto_old;

% memorizzo lo scarto nel vettore degli scarti

v_scarti(iter) = scarto;

% passo all'iterazione sucessiva

xold = xnew;

scarto_old = scarto;

end

R = -log10(asint);

end

ELIMINAZIONE DI GAUSS – MAIN

A = [1 3 5; 0 2 -3; 0 0 8];

b = [1; 2; 3];

[U,c] = gauss(A,b);

x = SosIndietro(U,c);

for i = size(b,1)

fprintf('x(%d): %f\n', i, x(i))

end

GAUSS – FUNCTION

% funzione per l'eliminazione di gauss

function [U,c] = gauss(A,b)

% U triangolare alta con diagonale unitaria

% c è il vettore risultato dopo l'eliminazione

[n,m] = size(A);

% verifica delle dimensioni della matrice

if n~=m

fprintf('errore: matrice non quadrata');

end

c = zeros(n,1);

U = eye(n);

for k = 1:n

for j = k+1:n

U(k,j) = A(k,j)/A(k,k);

end

c(k) = b(k)/A(k,k);

for i = k+1:n

for j = k+1:n

A(i,j) = U(k,j)*A(i,k)-A(i,j);

end

b(1) = c(k)*A(i,k)-b(i);

end

end

x(n) = c(n);

for i = n-1:-1:1

x(i) = c(i)-(U(i,(i+1):n)*x((i+1):n)');

end

end

SOSTITUZIONE INDIETRO – FUNCTION

function x = SosIndietro(U,c)

n = size(U,1);

x = zeros(n,1);

for i = n:-1:1

som = 0.0;

for j = i+1:n

som = som+U(i,j)*x(j);

end

x(i) = (c(i)-som)/U(i,i);

end

end

METODO SOR – MAIN

% apertura del file

fid = fopen('risul_sol.dat', 'w');

% chiamata alla funzione sor_comp

[xnew, iter, v_scarti, asint, R] = sor_comp(A, b, x0, omega, toll,

itmax);

% stampo

ST = [2:iter; v_scarti(2:iter)'; asint': R'];

fprintf(fid, '%5s %12s %12s\n', 'iter', 'scarto','asint', 'R');

fprintf(fid, '%5d, %15.6e, %15.6e, %15.6e\n', ST); fclose(fid);

% profilo di convergenza

semilogy(2:iter, abs(v_scarti(2:iter)), '-*b');

% calcolo il raggio spettrale di A

[~,n] = size(A);

Finv = diag(1./diag(A));

% calcolo triangolare bassa e alta

L = tril(A, -1);

U = triu(A, +1);

D = diag(diag(A));

E_omega = inv(D+omega*L)*((1-omega)*D-omega*U);

ragg_spett = max (abs(eig(E_omega)));

SOR_COMP - FUNCTION

function [xnew, iter, v_scarti, asint, R] = sor_comp(A, b, x0,

omega, toll, itmax)

% dimensione A

n = size (A, 1);

% preallocazione v_scarti, xnew

v_scarti = ones(itmax, 1);

xnew = zeros(n,1);

% inizializzo le variabili

xold = x0;

iter = 0;

scarto = 2.0*toll;

scarto_old = 1.0;

while scarto > toll && iter < itmax

iter = iter + 1;

for i = 1:n

% sommatoria parte bassa

som_L = 0.0;

for j = 1:i-1

som_L = som_L + A(i,j)*xnew(j);

end

% sommatoria parte alta

som_U = 0.0;

for j = i+1:n

som_U = som_U + A(i,j)*xold(j);

end

xseid(i) = (b(i)-som_L-som_U)/A(i,i);

xnew(i) = xold(i)+omega*(xseid(i)-xold(i));

end

scarto = norm(xnew-xold,2);

% memorizzo lo scarto nel vettore degli scarti

v_scarti(iter) = scarto;

% iterazione sucessiva

xold = xnew;

scarto_old = scarto;

end

% calcolo la costante di convergenza e la velocità

asint = v_scarti(2:iter)./v_scarti(1:iter-1);

R = -log10(asint);

end

SOR_MAT – FUNCTION

function [xnew, iter, v_scarti,