%基于LU分解的矩阵求逆
clc;
clear;
close all;
N = 6;
Source = randint(N,N,[0,10]);
Source_inv = zeros(N,N);
L = eye(N,N);
L_inv = eye(N,N);
U = zeros(N,N);
U_inv = zeros(N,N);
for k=1:1:N
if k==1
U(k,:) = Source(k,:);
else
for j=k:1:N
U(k,j) = Source(k,j) - L(k,1:(k-1))*U(1:(k-1),j);
end
end
if U(k,k)==0
error('Source error');
end
for i=k+1:1:N
L(i,k) = ( Source(i,k) - L(i,1:(k-1))*U(1:(k-1),k) )/U(k,k);
end
end
%inv(L)
for i=1:1:N
L_inv(i,i) = 1/L(i,i);
end
for i=2:1:N
for j=1:1:i-1
L_inv(i,j) = -L_inv(i,i)*(L(i,j:(i-1))*L_inv(j:(i-1),j));
end
end
%inv(U)
for i=1:1:N
U_inv(i,i) = 1/U(i,i);
end
for i=1:1:N-1
k=1;
for j=i+1:1:N
U_inv(k,j) = -U_inv(k,k)*(U(k,(k+1):j)*U_inv((k+1):j,j));
k=k+1;
end
end
Source_inv = U_inv*L_inv;
err = Source_inv - inv(Source);
err = err(:);
plot(err);