matlab_2_7

gaussSeidel.m

function [ x,iter ] = gauss_seidel(A,b,itermax,toll,x0)
%ESERCIZIO_2_7_F1 Summary of this function goes here
%   Detailed explanation goes here
% [ x,iter ] = gauss_seidel(A,b,itermax,toll,xo)
%input: 
%     A->matrice sistema lineare ,
%     b-> vettore termine noto,
%     itermax -> numero massimo di iterazioni,
%     toll-> tolleranza richiesta,
%     x0 -> vettore soluzione iniziale
%output:
%     x-> soluzione sistema lineare
%     iter -> numero iterazioni effettuate

%Ax=b -> (D+C)x=b -> Dx=b-Cx
D=tril(A); %gauss seidel: trangolare inferiore 
C=A-D;

n=length(b);
B= -inv(D)*A + eye(n); %B matrice di interazione
rho=max(abs(eig(B))); %raggio spettrale
disp('raggio spettrale')
disp(rho);

for iter=1:itermax
    x=D\(b-C*x0); %al primo passo uso x0
    if norm(x-x0,'inf') <= toll*norm(x,'inf') %condizione
        break; %esce dal ciclo for 
    end   
    x0=x;%poi x0 diventa la x di prima
end

%non devo fare nessun return, è utomatico e ritorna automaticamente i parametri descriti in output

end

jacobi.m

function [ x,iter ] = jacobi(A,b,itermax,toll,x0)
%ESERCIZIO_2_7_F1 Summary of this function goes here
%   Detailed explanation goes here
% [ x,iter ] = jacobi(A,b,itermax,toll,xo)
%input: 
%     A->matrice sistema lineare ,
%     b-> vettore termine noto,
%     itermax -> numero massimo di iterazioni,
%     toll-> tolleranza richiesta,
%     x0 -> vettore soluzione iniziale
%output:
%     x-> soluzione sistema lineare
%     iter -> numero iterazioni effettuate

%Ax=b -> (D+C)x=b -> Dx=b-Cx
D=diag(diag(A)); %jacobi: diagonale di A %diag di diag restituisce la diagonale di un vettore una matrice
C=A-D;

% x = invD * (b-Cx)
%   = invD*b -invD*C*x
%   = invD*b -invD*C*(invD*b - invD*C*x)
%   = ... + invD*C*invD*C*x
%   = ... + (-invD*C)^2*x
%   = ... + (-invD*C)^3*x
%   = ... + B^3*x
% B = -invD*C = -invD*(A-D) = -invD*A +invD*D 
%   = -invD*A + I
n=length(b);
B= -inv(D)*A + eye(n); %B matrice di interazione
rho=max(abs(eig(B))); %raggio spettrale
disp('raggio spettrale: ')
disp(rho);

for iter=1:itermax
    x=D\(b-C*x0); %al primo passo uso x0
    if norm(x-x0,'inf') <= toll*norm(x,'inf') %condizione
        break; %esce dal ciclo for 
    end   
    x0=x;%poi x0 diventa la x di prima
end

%non devo fare nessun return, è utomatico e ritorna automaticamente i parametri descriti in output

end

matlab.m

%Esercizio_2_7
%vedere i file funzioni

clear all
clc

itermax=100;
toll=1e-7; %1*10^-7

A=[1 -2 2; -1 1 -1; -2 -2 1];
b=sum(A,2);

n=length(b); %dimensione del problema
x0=zeros(n,1);

disp('Jacobi')
[xJ,iterJ] = jacobi(A,b,itermax,toll,x0)
disp('Gauss Seidel')
[xGS,iterGS] = gauss_seidel(A,b,itermax,toll,x0)


%% Esercizio con matrice A3

A3=[2 -1 1; 2 2 2; -1 -1 2];
b3=sum(A3,2);

n=length(b); %dimensione del problema

disp('Jacobi')
[xJ,iterJ] = jacobi(A3,b3,itermax,toll,x0)
disp('Gauss Seidel')
[xGS,iterGS] = gauss_seidel(A3,b3,itermax,toll,x0)


%% Esercizio con matrice A2

A2=[4 0 2/5;0 5 2/5; 5/2 2 1];


%% Esercizio con matrice A4