矩阵的三角分解算法实验报告

作者：戴宇童

如有任何问题请联系 [e-mail]

LR 分解,LDR 分解算法

LR分解的代码

function [A,L,R]=myLR(A,op)%op是一个数值，op=1表示执行节约内存版的计算方式，否则执行普通LR分解%如果执行op=1，则所求的单位下三角存在A的下三角部分（不含对角线），所求的上三角矩阵R存储在A的上三角部分（包含对角线）n=size(A,2);for i=1:n    Z=det(A(1:i-1,1:i-1));    if Z==0         fprintf('Factorization is impossible')    endendif op==1    A(1,:)=A(1,:);    A(2:end,1)=A(2:end,1)/A(1,1);    for r=2:n        for j=r:n        A(r,j)=A(r,j)-A(r,1:r-1)*A(1:r-1,j);        end        for i=r+1:n        A(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r))/A(r,r);        end    endelse    R=zeros(size(A));    L=eye(size(A));    R(1,:)=A(1,:);    L(2:end,1)=A(2:end,1)/A(1,1);for r=2:n    for j=r:n        R(r,j)=A(r,j)-L(r,1:r-1)*R(1:r-1,j);    end    for i=r+1:n        L(i,r)=(A(i,r)-L(i,1:r-1)*R(1:r-1,r))/R(r,r);    endendend

LDR分解代码

function [L,D,Rn]=myLDR(A,op)%op是一个数值，op=1表示执行节约内存版的计算方式，否则执行普通分解%{如果执行节约内存的计算方式，那么得到的矩阵仍然记作A。A的下三角（不含对角线）是L的下三角元素（不含对角线）A的对角元是对角阵D的对角元A的上三角（不含对角线）是R的上三角元素（不含对角线）%}n=size(A,2);if op==1    A=myLR(A,1);    for i=1:n-1        for j=i+1:n            A(i,j)=A(i,j)/A(i,i);        end    end    Aelse    [A,L,R]=myLR(A,2);    for i=1:n    D(i,i)=R(i,i);    end    for i=1:n    Rn(i,i)=1;    end    for i=1:n-1        for j=(i+1):n        Rn(i,j)=R(i,j)/R(i,i);        end    endend

例子

对A1=⎛⎝⎜⎜21131249−6⎞⎠⎟⎟进行LR分解与LDR分解

利用代码

A=[2,3,4;1,1,9;1,2,-6];myLR(A,1)

得到结果为

利用代码

A=[2,3,4;1,1,9;1,2,-6];[A,L,R]=myLR(A,2)

得到结果为

利用代码

A=[2,3,4;1,1,9;1,2,-6];myLDR(A,1)

得到结果为

利用代码

A=[2,3,4;1,1,9;1,2,-6];[L,D,Rn]=myLDR(A,2)

得到结果为

Cholesky 分解

算法一

%Cholesky分解(书上算法一)%函数使用说明A是输入待分解的矩阵%op=1表示执行节约内存版的计算方式%如果使用op=1,则所求的下三角矩阵存储在A的下三角部分（包括对角线）function [A,T]=myCHOL1(A,op)if det(A)==0;    fprint('Fractorization is impossible')endn=size(A,2);if op==1    A(1,1)=sqrt(A(1,1));    for i=2:n        A(i,1)=A(1,i)/A(1,1);    end    for i=2:n        s=0;        for k=1:i-1            s=s+A(i,k)^2;        end        A(i,i)=sqrt(A(i,i)-s);        for j=i+1:n            t=0;            for k=1:i-1                t=t+A(i,k)*A(j,k);            end            A(j,i)=(A(i,j)-t)/A(i,i);        end    endelse    T(1,1)=sqrt(A(1,1));    for j=2:n        T(j,1)=A(1,j)/sqrt(A(1,1));    end    for i=2:n        sum=0;        for k=1:(i-1)            sum=sum+T(i,k)^2;        end        T(i,i)=sqrt(A(i,i)-sum);        for j=i+1:n            sum=0;            for k=1:(i-1)            sum=sum+T(i,k)*T(j,k);            end        T(j,i)=(A(i,j)-sum)/T(i,i);        end        for j=1:(i-1)            T(j,i)=0;        end    endend

算法二

%Cholesky分解(书上算法二)%函数使用说明A是输入待分解的矩阵%输出矩阵T是下三角正交阵%A=[4,6,10;6,58,29;10,29,38]function myCHOL3(A)if det(A)==0 || det(A)<0;    fprintf('Fractorization is impossible')endn=size(A,2);for k=1:n    for j=k:n        T(j,k)=A(k,j)/sqrt(A(k,k));    end        for i=k+1:n            for t=k+1:n            A(i,t)=A(i,t)-A(i,k)*A(k,t)/A(k,k);            end        endendT

例子

对A2=⎛⎝⎜⎜120222023⎞⎠⎟⎟进行Cholesky分解

算法一

利用代码

A2=[1,2,0;2,2,2;0,2,3];myCHOL1(A2,1)

得到结果
矩阵A2不可分解，究其原因是因为它不正定！

为了验证算法的正确性，我们选取一个正定阵来验证

A3=⎛⎝⎜⎜2003491451413⎞⎠⎟⎟

利用代码

A3=[4,6,10;6,58,29;10,29,38];myCHOL1(A3,1)

得到结果

利用代码

A3=[4,6,10;6,58,29;10,29,38];[A3,T]=myCHOL1(A3,2)

得到结果

算法二

利用代码

A3=[4,6,10;6,58,29;10,29,38];myCHOL3(A3)

得到结果

求解线性方程组Ax=b

要求A是方阵且可以作LR分解

%求解线性方程组Ax=bfunction x=mysolve(A,b)n=size(A,2);[A,L,R]=myLR(A,2);%求解Lb*=bb_star(1)=b(1)/L(1,1);for i=2:n    s=0;    for k=1:i-1        s=s+b_star(k)*L(i,k);    end    b_star(i)=(b(i)-s)/L(i,i);end%求解Rx=b_starx(n)=b_star(n)/R(n,n);for i=1:n-1    t=0;    for k=n-i+1:n        t=t+x(k)*R(n-i,k);    end    x(n-i)=(b_star(n-i)-t)/R(n-i,n-i);end

例子

设:
A=[2,0,1;0,1,0;1,0,2];
b=[3;-4;0];

利用命令

 x=mysolve(A,b)

可以得到答案

求解上三角或下三角矩阵的逆

%求解上三角矩阵R或者下三角矩阵L的逆%算法:L(b1,...,bn)=(e1,...,en) %R(r1,...,rn)=(e1,...,en)function INV=myinverse(A)I=eye(size(A));n=size(A,2);if A(1,2)==0        INV=zeros(size(A));    for j=1:n    INV(1,j)=I(1,j)/A(1,1);    for i=2:n        s=0;        for k=1:i-1        s=s+INV(k,j)*A(i,k);        end    INV(i,j)=(I(i,j)-s)/A(i,i);    end    endelse        for j=1:n    INV(n,j)=I(n,j)/A(n,n);        for i=1:n-1            t=0;            for k=n-i+1:n                t=t+INV(k,j)*A(n-i,k);            end        INV(n-i,j)=(I(n-i,j)-t)/A(n-i,n-i);        end    endend

例子

假设:
A1=[1,0;2,3];
A2=[1,2;0,3]:

利用命令

myinverse(A1)myinverse(A2)

可以得到结果