Esercizi e prove svolte della parte Matlab di Informatica

INTERPOLAZIONE

Interpol.mx = linspace(0,20,100);y = cos(x);p = polyfit(x,y,3);yp = polyval(p,x);plot(x,y,'x');hold on;plot(x,yp);

Interpolreale.mx = linspace(0,2,100);y = cos(x);p = polyfit(x,y,4);yp = polyval(p,x);plot(x,y,'x');hold on;plot(x,yp);

Interpolesatta.mx = linspace(0,20,10);y = cos(x);p = polyfit(x,y,9);xp = linspace(0,20,1000);yp = polyval(p,xp);plot(x,y,'x');hold on;plot(xp,yp);

TRAPEZI

TrapSimple.m

function It = TrapSimple(a,b,f)
    fa = f(a);
    fb = f(b);
    It = (b-a)*(fa+fb)/2;
end

TrapComp.m

function It = TrapComp(a,b,n,f)
    h = (b-a)/n;
    x_sx = a;
    It = 0.0;
    for i = 1:n
        x_dx = x_sx + h;
        It_tmp = TrapSimple(x_sx,x_dx,f);
        It = It + It_tmp;
        x_sx = x_dx;
    end
end

TrapComp_2.m

function It = TrapComp_2(a,b,n,f)
    x_i = linspace(a,b,n+1);
    f_i = f(x_i);
    h = (b-a)/n;
    It = h*(0.5*(f_i(1) + f_i(end)) + sum(f_i(2:end-1)));
end

main.m

clear
clc

fun = @(x) sin(x);
a = 0.0;
b = pi/2;
n = 10;

Itsimple = TrapSimple(a,b,fun);
fprintf('Integrale semplice = %e\n',Itsimple);

Itcomp = ...

TrapComp(a,b,n,fun);
fprintf('Integrale composto 1 = %e\n',Itcomp);
Itcomp_2 = TrapComp_2(a,b,n,fun);
fprintf('Integrale composto 2 = %e\n',Itcomp_2);
Iv = -cos(b) + cos(a);
fprintf('Integrale vero = %e\n',Iv);

maintempo.mclearclcfun = @(x) sin(x);
a = 0.0;
b = pi/2;
n = 10;
Itsimple = TrapSimple(a,b,fun);
fprintf('Integrale semplice = %e\n',Itsimple);

tic;
Itcomp = TrapComp(a,b,n,fun);
toc
fprintf('Integrale composto 1 = %e\n',Itcomp);

tic;
Itcomp_2 = TrapComp_2(a,b,n,fun);
toc
fprintf('Integrale composto 2 = %e\n',Itcomp_2);

Iv = -cos(b) + cos(a);
fprintf('Integrale vero = %e\n',Iv);

mainerrore.mclearclcfun = @(x) sin(x);
a = 0.0;
b = pi/2;
k = 10;
n = 1;
Iv = -cos(b) + cos(a);
for i = 1:k
    Itcomp_2 = TrapComp_2(a,b,n,fun);
    Et = abs(Itcomp_2 - Iv);
    fprintf('n: %d It %15.6e Et %15.6e\n',n,Itcomp_2,Et);
    n = n*2;
end

mainerrorerapporto.mclearclcfun = @(x) sin(x);
a = 0.0;
b = pi/2;
k = 10;
Iv = -cos(b) + cos(a);
n = 1;
Et_old = 1.0;
for i = 

ratio%10.2f\n',n,Qcomp,Ir,Eq,Eq_old/Eq);n = n*2;Eq_old = Eq;Qcomp_old = Qcomp;end

CavSimpComp.m
function Iq = CavSimpComp(a,b,n,f)
    x_i = linspace(a,b,2*n+1);
    f_i = f(x_i);
    h = (b-a)/n;
    pari = 2:2:2*n;
    dispari = 3:2:2*n-1;
    Iq = (h/6)*((f_i(1) + f_i(end)) + 4*sum(f_i(pari)) + 2*sum(f_i(dispari)));
end

METODI DIRETTI ITERATIVI STAZIONARI – BACKFORWSOLVEmain.m
clear
clc

% 1° metodo di soluzione
A = [5 3 1;1 4 2;3 1 6];
b = [9; 7; 10];
x1 = A\b;
fprintf('Soluzione 1: \n');
fprintf('%10.2f\n',x1);

% 2° metodo di soluzione
[L,U] = lu(A);
% Risolvo Ly = b
y = L\b;
% Risolvo Ux = y
x2 = U\y;
fprintf('Soluzione 2: \n');
fprintf('%10.2f\n',x2);

% Calcolo soluzione con BackForwSolve fatta da mex
x3 = BackForwSolve(L,U,b);
fprintf('Soluzione 3: \n');
fprintf('%10.2f\n',x3);

BackForwSolve.m
function x = BackForwSolve(L,U,b)
    n = size(L,1);
    y = zeros(n,1);
    x = zeros(n,1);
    
    % Sostituzione avanti con L
    for i = 1:n
        som = 0.0;
        for j = 1:i-1
            som = som + L(i,j)*y(j);
        end
        y(i) = (b(i) - som)/L(i,i);
    end
    
    % Sostituzione indietro con U
    for i = n:-1:1
        som = 0.0;
        for j = i+1:n
            som = som + U(i,j)*x(j);
        end
        x(i) = (y(i) - som)/U(i,i);
    end
end

L(i,j)*y(j);
end
y(i) = b(i) - som;
end

% Sostituzione indietro con U
for i = n:-1:1
som = 0.0;
for j = i+1:n
som = som + U(i,j)*x(j);
end
x(i) = (y(i) - som)/U(i,i);
end
end

JACOBI

Jacobi.m
clear
clc

A = [9 1 2;2 8 1;-1 1 6];
b = [12;11;6];
itmax = 100;
toll = 1.e-5;
x0 = zeros(3,1);
n = size(A,1);

% Inizio Jacobi
iter = 0;
scarto = 2 * toll;
scarto_old = scarto;
x_old = x0;
x_new = zeros(n,1);

while scarto > toll && iter < itmax
iter = iter + 1;
for i = 1:n
som1 = 0;
for j = 1:i-1
som1 = som1 + A(i,j) * x_old(j);
end
som2 = 0;
for j = i+1:n
som2 = som2 + A(i,j) * x_old(j);
end
x_new(i) = (b(i) - som1 - som2) / A(i,i);
sc = x_new - x_old;
scarto = norm(sc);
end
rad = scarto / scarto_old;
x_old = x_new;
scarto_old = scarto;
fprintf('iter: %d | scarto %e | raggio %e\n',iter,scarto,rad);
end

Jacobi_2.m
clear
clc

A = [9 1 2;2 8 1;-1 1 6];
b = [12;11;6];
itmax = 100;
toll = 1.e-5;
x0 = zeros(3,1);
n = size(A,1);

% Inizio Jacobi_2
iter = 0;
scarto = 2 * toll;
scarto_old = scarto;
x_old = x0;
x_new = zeros(n,1);

% xkp1 = xk +

D^-1*rk = Ej*xk + q% con Ej = I - D^-1*A e q = D^-1*b

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

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

while scarto > toll && iter < itmax

iter = iter + 1;

x_new = Ej * x_old + q;

for i = 1:n

som1 = 0;

for j = 1:i-1

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

end

som2 = 0;

for j = i+1:n

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

end

x_new(i) = (b(i) - som1 - som2) / A(i,i);

sc = x_new - x_old;

scarto = norm(sc);

end

scarto = norm(x_new - x_old);

rad = scarto/scarto_old;

x_old = x_new;

scarto_old = scarto;

fprintf('iter: %d | scarto %e | raggio %e\n',iter,scarto,rad);

end

GAUSS-SEIDELSeidel.m

clear

clc

A = [9 1 2;2 8 1;-1 1 6];

b = [12;11;6];

itmax = 100;

toll = 1.e-5;

x0 = zeros(3,1);

n = size(A,1);

% Inizio Seidel

iter = 0;

scarto = 2 * toll;

scarto_old = scarto;

x_old = x0;

x_new = zeros(n,1);

while scarto > toll && iter < itmax

iter = iter + 1;

for i = 1:n

som1 = 0;

for j = 1:i-1

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

end

som2 = 0;

for j = i+1:n

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

end

x_new(i) = (b(i) - som1 - som2) / A(i,i);

sc = x_new - x_old;

scarto = norm(sc);

end

scarto = norm(x_new - x_old);

rad = scarto/scarto_old;

x_old = x_new;

scarto_old = scarto;

fprintf('iter: %d | scarto %e | raggio %e\n',iter,scarto,rad);

end

A(i,i);sc = x_new - x_old;scarto = norm(sc);endrad = scarto/scarto_old;x_old = x_new;scarto_old = scarto;fprintf('iter: %d | scarto %e | raggio %e\n',iter,scarto,rad);end Seidel_2.m

clear

clc

A = [9 1 2;2 8 1;-1 1 6];

b = [12;11;6];

itmax = 100;

toll = 1.e-5;

x0 = zeros(3,1);

n = size(A,1);

% Inizio Seidel_2

iter = 0;

scarto = 2 * toll;

scarto_old = scarto;

x_old = x0;

x_new = zeros(n,1);

% xkp1 = xk + (L+D)^-1*rk = Es*xk + q

% con Es = I - (L+D)^-1*A e q = (L+D)^-1*b

Es = eye(n) - inv( tril(A) )*A;

q = inv( tril(A) )*b;

while scarto > toll && iter < itmax

iter = iter + 1;

xnew = Es*x_old + q;

scarto = norm(xnew - x_old);

rad = scarto/scarto_old;

x_old = x_new;

scarto_old = scarto;

fprintf('iter: %d | scarto %e | raggio %e\n',iter,scarto,rad);

end

ESERCIZI TEORICI

NORMA EUCLIDEA

NormaEuc.m

clear

clc

n = 10;

E = rand(n);

r = rand(n,1);

% Metodo 1

% s = (I + Et * E)*r = r + Et * v con v = E*r

v1 = E*r;

s1 = r + E'*v1;

norma_1 = norm(s1);

% Metodo 2

v2 = matvet(E,r);

s2 = r +

matvet(E',v2);
norma_2 = neuc(s2);
fprintf('Ris 1: %e\n',norma_1);
fprintf('Ris 2: %e\n',norma_2);
matvet.m
function y = matvet(A,x)
[n,m] = size(A);
y = zeros(n,1);
for i = 1:n
    for j = 1:m
        y(i) = y(i) + A(i,j)*x(j);
    end
end
end
neuc.m
function alpha = neuc(x)
n = size(x,1);
alpha = 0;
for i = 1:n
    alpha = alpha + x(i)^2;
end
alpha = sqrt(alpha);
end
NORMA ASSOLUTA DI UN VETTORE
dati.txt
3
e3
n10
1
2
B
1
x
-1
3
1
1
-1
2
2
1
1.e-4
toll
100
itmax
main.m
clear
ifile = fopen('dati.txt','r');
l = fscanf(ifile,'%e',1); % Legge la potenza
n = fscanf(ifile,'%e',1); % Legge la dimensione
B = fscanf(ifile,'%e',[n,n]); % Legge la matrice B
B = B';
x = fscanf(ifile,'%e',n); % Legge x
for i = 1 : l
    z = matvet(B,x);
    x = z;
end
assoluta = normass(x);
fprintf('Norma assoluta di x: %e\n',assoluta);
matvet.m
function v = matvet(A,x)
nr = size(A,1);
nc = size(A,2);
for i = 1 : nr
    v(i,1) = 0.0;
    for j = 1 : nc
        v(i,1) = v(i,1) + A(i,j) * x(j,1);
    end
end
end

normass.m

function norma = normass(x)
    n = size(x,1);
    norma = 0;
    for i = 1 : n
        norma = norma + abs(x(i,1));
    end
end

METODO DELLE POTENZE PER AUTOVALORE MAX

dati.txt
3
n
10
1
2
A
1
x0
-1
3
1
1
-1
2
2
1
1.e-4
toll
100
itmax
main.m
clear

ifile = fopen('dati.txt','r');
n = fscanf(ifile,'%e',1); % Legge la dim
A = fscanf(ifile,'%e',[n,n]); % Legge matr A
A = A';
z0 = fscanf(ifile,'%e',n); % Legge z0, toll, itmax
toll = fscanf(ifile,'%e',1);
itmax = fscanf(ifile,'%e',1);
iter = 0;
residuo = 2 * toll;
z_old = z0;
while residuo > toll && iter < itmax
    iter = iter + 1;
    % Tre passi delle potenze
    y = matvet(A,z_old);
    lambda = eucl(y);
    z_new = y / lambda;
    % Calcolo residuo
    num = eucl(y - lambda * z_old);
    den = eucl(z_old);
    residuo = num / den;
    % Passo all'iterazione successiva
    z_old = z_new;