Resolución de Sistemas Lineales en MATLAB: Factorización LU, Cholesky y Análisis de Errores

Resolución de sistemas lineales en MATLAB

A continuación se presentan varios scripts y funciones en MATLAB para la resolución de sistemas de ecuaciones lineales mediante factorización LU, descomposición de Cholesky y métodos de sustitución, junto con el cálculo de residuos, errores relativos y números de condición.

1. Sistema lineal con factorización LU (sin pivoteo)

clear all
format short e

a=ones(10)+eye(10);
a(1,10)=0.5;

for i=1:10
    b(i,1)=2*i;
end

[l,u,ifallo]=gauslu(a);

if ifallo==0
    [y,ifallo]=solvetrinf(l,b);
    if ifallo==0
        [x,ifallo]=solvetrsup(u,y);
    end
end

if ifallo==0
    fprintf('Norma del residuo %14.6e \n',norm(a*x-b))
    fprintf('Solución calculada %14.6e \n',x(1))
    for i=2:10
        fprintf('%34.6e \n',x(i))
    end
end

2. Sistema con factorización LU con pivoteo por filas

clear all
format short e

n=10;
c=ones(n);

for j=1:n
    c(n+1-j,j)=2;
end
c(n,n)=0.5;

for j=1:n
    d(j,1)=2*(n+1-j);
end

% Factorización con vector de permutación
[l,u,p]=lu(c,'vector');

[y,ifallo]=solvetrinf(l,d(p));

if ifallo==0
    [x,ifallo]=solvetrsup(u,y);
end

if ifallo==0
    fprintf('Residuo %14.6e \n',norm(c*x-d))
    fprintf('Error relativo %14.6e \n',norm(xexacta-x)/norm(xexacta))
else
    fprintf('No éxito')
end

3. Construcción de la matriz A

A=diag(4*ones(n,1))+diag(ones(n-1,1),1)+diag(ones(n-1,1),-1);

4. Sistema tridiagonal con análisis de condicionamiento

a=ones(n-1,1);
b=4*ones(n,1);
c=a;
d=(1:n)';

A=...   % Matriz tridiagonal asociada

tiempoi=cputime;
[x_tri,ifail]=tridisolve2(a,b,c,d);
tiempo=cputime-tiempoi;

whos;

% Cotas superior e inferior del error relativo
COTSUP=cond(A)*norm(A*x_tri-d)/norm(d);
COTINF=(1/cond(A))*norm(A*x_tri-d)/norm(d);

5. Resolución mediante LU con medición de tiempo

A=...
d=...

tiempoi=cputime;
[l,u,p]=lu(A,'vector');
[y,ifallo]=solvetrinf(l,d(p));

if ifallo==0
    [x,ifallo]=solvetrsup(u,y);
    tiempo=cputime-tiempoi;
end

whos
COTSUP
COTINF

6. Construcción de una matriz particular

A=tril(-ones(n));
for i=1:n
    A(i,i)=1;
end
for i=1:n
    A(i,n)=1;
end

7. Función para resolver sistemas de Hilbert mediante Cholesky

function [error, residuo, ifallo]=hilbertsolve(n)

error=[];
residuo=[];

A=hilb(n);
xexacta=ones(n,1);
b=A*xexacta;

% Factorización de Cholesky: A = L*L'
[L,ifallo]=chol(A,'lower');
Ltras=L';

if ifallo==0
    [y,ifallo]=solvetrinf(L,b);
    if ifallo==0
        [x,ifallo]=solvetrsup(Ltras,y);
        x
    else
        disp('Error')
        return
    end
else
    disp('Error')
    return
end

% Cálculo del error relativo y del residuo
error=norm(x-xexacta)/norm(x);
residuo=norm(b-A*x);

end

8. Cálculo de la inversa mediante LU

% Resolver A*C = D, aplicando LU a cada columna de la identidad
d=eye(n);

for k=1:n
    b=d(:,k);
    [l,u]=lu(a);
    if ifallo==0
        [y,ifallo]=solvetrinf(l,b);
        if ifallo==0
            [x,ifallo]=solvetrsup(u,y);
        end
    end
    c(:,k)=x;
end