在稀疏矩阵网站上找到了一个稀疏矩阵，准备对比高斯消去，LU分解和Cholesky分解的运行效率，但是一直报错，孩子改不出来了，呜呜呜呜呜！求大神解答
LU分解法
load('bcspwr03.mat')
A=full(Problem.A);
b=rand(118,1);
A
function [L,U]=myLU(A)
[n,n]=size(A);
L=zeros(n,n);
U=zeros(n,n);
%% 对A的LU分解
for i=1:n
L(i,i)=1;
end
for k=1:n
for j=k:n
U(k,j)=A(k,j)- sum(L(k,1:k-1).*U(1:k-1,j)');
end
for i=k+1:n
L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k);
end
end
%% 求解y过程
%消去过程
function x=GaussianSolver(L,b)
[n,k] = size(A);
[n,k] = size(L);
x=zeros(n,1);
for j = 1:n-1
for i = j+1:n
mul = L(i,j)/L(j,j);
L(i,:)= L(i,:)-mul*L(j,:);
b(i)= b(i)-mul*b(j);
end
end
%回代过程
for i=n:-1:1
sum=0;
for j=n:-1:i+1
sum=sum+y(j)*L(i,j);
end
y(i)=(b(i)-sum)/L(i,i);
end
end
%% 求解x过程
%消去过程
for j = 1:n-1
for i = j+1:n
mul = U(i,j)/U(j,j);
U(i,:)= U(i,:)-mul*U(j,:);
y(i)= y(i)-mul*y(j);
end
end
%回代过程
for i=n:-1:1
sum=0;
for j=n:-1:i+1
sum=sum+x(j)*U(i,j);
end
x(i)=(y(i)-sum)/U(i,i);
end
end
Cholesky分解法
load('bcspwr03.mat')
A=full(Problem.A)
b=rand(118,1);
function [A,b]=creatMaxtrix(n)
T=zeros(n-1,n-1);
for i=1:n-1
T(i,i)=2;
end
for i=1:n-2
T(i,i+1)=-1;
T(i+1,i)=-1;
end
I=eye(n-1,n-1);
k=(n-1)^2;
A=zeros(k,k);
b=ones(k,1);
j=1;
for i=1:n-1
A(j:j+n-2,j:j+n-2)=T+2*I;
j=j+n-1;
end
j=1;
for i=1:n-2
A(j:j+n-2,j+n-1:j+2*n-3)=-I;
A(j+n-1:j+2*n-3,j:j+n-2)=-I;
j=j+n-1;
end
end
function []=Cholesky(n)
[A,b]=creatMaxtrix(n);
[n,n]=size(b);
x=A\b;
%% Cholesky分解
for k=1:n
A(k,k)=sqrt(A(k,k));
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
for j=k+1:n
A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);
end
end
%% 求解
%解下三角形方程组：前代法
L=tril(A,0);
for j=1:n-1
b(j)=b(j)/L(j,j);
b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);
end
b(n)=b(n)/L(n,n);
%解上三角形方程组：回代法
U=L';
for j=n:-1:2
b(j)=b(j)/U(j,j);
b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);
end
b(1)=b(1)/U(1,1);
%用二范数来衡量误差
norm(x-b)
end
求求了
不知道你这个问题是否已经解决, 如果还没有解决的话: