中科院矩阵分析与应用大作业

中科院矩阵分析与应用大作业

按照上面建立函数文件，脚本文件。名字不要随便改，要改的话要同内部函数名一起修改。下面放代码，总共七段代码，第一个是主程序进行人机交互，其余为接口函数实现矩阵的各个分解功能。另外容我小小的收点下载券，毕竟花了大概8小时左右。

%%矩阵分解：LU分解；QR分解(schmidt)；Householder；Givens %%2016.11.22

clc,clear all;

Matrix=input('请输入需要进行分解的矩阵：');

disp('请选择需要对矩阵进行的操作:');

disp('1.LU分解');disp('2.QR分解');

disp('3.Householder约简');disp('4.Givens约简');

n=input('');

clc;%清除命令窗口内容

Matrix%显示输入矩阵

%%分解操作

switch n

case1

LUfactor(Matrix);%LU分解

case2

QR(Matrix);%QR分解

case3

houseHolder(Matrix);%householder约简case4

Givens(Matrix);%Givens约简

end

function LUfactor(A)

[m n]=size(A);

%判断矩阵是否为方阵

if m~=n

error('矩阵非方阵，请重新输入！'); end

%判断矩阵是否奇异

if det(A)==0

error('矩阵奇异，请重新输入！');

end

%判断矩阵顺序主子式是否全不为0

n=size(A,1);

flagMat=zeros(n,1);

for i=1:n

if det(A(1:i,1:i))==0

flagMat(i)=1;

end

end

%以顺序主子式是否含0来决定采用的LU分解方式if any(flagMat)==0

disp('顺序主子式均不为0，采用A=LU');

disp('***LU分解结果***')

[L U]=LUFull(A)

else

disp('顺序主子式存在0，采用PA=LU');

disp('***LU分解结果***')

[L U P]=LUPart(A)

end

%对矩阵A进行LU分解(完全主元法)

%L,U矩阵为输出变量；A矩阵为输入变量

U=zeros(size(A));%对U尽可能初始化

U(1,:)=A(1,:);

L=eye(size(A));%对L尽可能初始化

L(2:end,1)=A(2:end,1)/A(1,1);

n=size(A,1);

for i=2:n%Dolittle公式

for j=i:n

U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);

end

for k=i:n

L(k,i)=(A(k,i)-L(k,1:i-1)*U(1:i-1,i))/U(i,i);

end

end

%对矩阵A进行LU分解(部分主元法)

%L,U,P矩阵为输出变量，A矩阵为输入变量

%对输入矩阵A重构

n=size(A,1);

B=zeros(n,1);

for i=1:n

B(i,1)=B(i,1)+i;

end

A=[A B];

%目标矩阵

tempMat=zeros(size(A));

for j=1:n%第j列

[~,row]=max(abs(A(j:end,j)));%找出第j列最大元素,返回所在行A([row+j-1j],:)=A([j row+j-1],:);%交换最大元素行与第j行

tempMat(j,:)=A(j,:);%将第j行元素复制到tempMat第j行for i=j:n%第i行

if i+1<=n

tempMat(i+1,j)=A(i+1,j)/A(j,j);

end

end

%生成新的A(高斯消去)

for k=j:n%第k行

if k+1<=n

ratio=-A(k+1,j)/A(j,j);

A(k+1,j+1:n)=ratio*A(j,j+1:n)+A(k+1,j+1:n);

end

end

end

%提取U矩阵

U=triu(tempMat(:,1:n),0);%upper triangle

%提取L矩阵

L=tril(tempMat(:,1:n),0);%lower triangle

for i=2:n

for j=1:(i-1)

L(i,j)=L(i,j)./L(j,j);

end

end

%完善L矩阵的对角线

for i=1:n

L(i,i)=1;

end

%提取P矩阵

tempP=tempMat(:,n+1); P=zeros(n,n);

for i=1:n

P(i,tempP(i))=1; end

function QR(A)

%%QR分解之施密特正交

%%2016.11.21

%%判断是否能进行QR分解

[rows,cols]=size(A);

if rank(A)~=cols

error('该矩阵无法进行QR分解！');%检测参数异常,停止程序end

%%变量预分配类存

Q=zeros(rows,cols);

R=zeros(cols,cols);

uk=zeros(rows,1);

%%施密特正交求Q矩阵

for i=1:1:cols

if i==1

Q(:,1)=A(:,1)/sqrt(A(:,1)'*A(:,1));

else

j=i-1;

while j>0

uk=Q(:,j)'*A(:,i)*Q(:,j)+uk;

j=j-1;

end

Q(:,i)=(A(:,i)-uk)/sqrt((A(:,i)-uk)'*(A(:,i)-uk));

uk=0;

end

end

%%求R矩阵

for i=1:1:cols

for j=i:1:cols

R(i,j)=Q(:,i)'*A(:,j);

end

end

disp('***schmidt正交QR分解结果***')

Q,R