%%%程序编写者 西北工业大学自动化学院 Email: yincwxa2013@mail.nwpu.edu.cn
%% All rights reserved
clear
clc
a=input('输入对称正定矩阵a=')
b=input('输入自由项b=')
n=length(a(:,1));
for k=1:n
if (det(a(1:k,1:k))<=0)
input('矩阵不是正定矩阵,请重新运行程序')
end
end
%%%%%%%%%%%分解A=L*L'
for i=1:n
t=0;
for s=1:i-1
t=t+L(i,s)^2;
end
L(i,i)=sqrt(a(i,i)-t);
for k=i+1:n
tt=0;
for s=1:i-1
tt=tt+L(i,s)*L(k,s);
end
L(k,i)=(a(k,i)-tt)/L(i,i);
end
end
%%%%%%%%%%%%%%%分解AX=b为Ly=b Lx=y
%%%%%%%%%%求y
for i=1:n
ttt=0;
for k=1:i-1
ttt=ttt+L(i,k)*y(k);
end
y(i)=(b(i)-ttt)/L(i,i);
end
%%%%%%%%%%求x
for i=n:-1:1
tttt=0;
for k=i+1:n
tttt=tttt+L(k,i)*x(k);
end
x(i)=(y(i)-tttt)/L(i,i);
end
x
转载本文请联系原作者获取授权,同时请注明本文来自殷春武科学网博客。
收藏
分享
分享到: