数值分析算法 MATLAB 实践 线性方程组 分解法

Lu分解法

% LUmethod分解矩阵
function [L,U]=LUmethod(A)
[rows,~]=size(A);
temp_mat=A;
L=zeros(rows);
for i=1:rows
    coefficient=temp_mat(:,i);
    coefficient=coefficient./coefficient(i);
    coefficient(1:i)=0;
    L(:,i)=coefficient;
    temp_mat=-coefficient*temp_mat(i,:)+temp_mat;
end
U=temp_mat;
L(eye(rows)==1)=1;
end
% LUmethod分解矩阵 求解线性方程组
function x=LUsolve(L,U,b)
[rows,~]=size(L);
aug_mat=[L,b];
for i=1:rows
    aug_mat(i,:)=aug_mat(i,:)./aug_mat(i,i);
    coefficient=-aug_mat(:,i);
    coefficient(1:i)=0;
    aug_mat=coefficient*aug_mat(i,:)+aug_mat;
end
aug_mat=[U,aug_mat(:,rows+1:end)];
for i=rows:-1:1
    aug_mat(i,:)=aug_mat(i,:)./aug_mat(i,i);
    coefficient=-aug_mat(:,i);
    coefficient(i:end)=0;
    aug_mat=coefficient*aug_mat(i,:)+aug_mat;
end
x=aug_mat(:,rows+1:end);
x=x';
end

function solution=LuFunmethon(M, Presion) 
 % LU分解 M为用户输入的增广矩阵
 % Precision为用户所输入的精度要求
    if nargin==2
      try
         digits(Precision);
      cath
         disp('你输入的精度有误');digits(10);
      end
    else
      digits(10);
    end     

    A=vpa(M)
    row=size(A,1);
    col=size(A,2);
    if ndims(A)~=2|(col-row)~=1
        disp('矩阵的大小有误');
        return
    end
    if det(M(:,1:row))==0
        disp('该方程的系数矩阵行列式为零');
        return
    end
%% 调用系统的LU命令
       [L,U,P]=lu(double(A));
      %% 回代求解过程
    for i=row:-1:1
        temp=U(i,col);
        for k=i+1:row
            temp=vpa(temp-t_solution(k)*U(i,k));
        end
            t_solution(i)=vpa(temp/U(i,i));
    end
    for i=1:row
        temp=t_solution(i);
        for k=1: i-1
            temp=vpa(temp-t_solution(k)*U(i,k));
        end
             solution(i)=temp;
    end
end

Cholesky分解

%% 0.平方根法解线性方程组，输出L矩阵和根
%% 1.对称正定矩阵的Cholesky分解
%对称正定矩阵A存在唯一的对角元素均为正数的下三角矩阵L，使得A=L*L'
%这种分解叫做Cholesky分解
A=[3,3,5;3,5,9;5,9,17];
b=[0;-2;-4];
%L=chol(A,'lower')基于矩阵A的对角线和下三角形生成下三角矩阵L，满足方程L*L'=A
L=chol(A,'lower')
%% 2.由Ly=b得到y
y=L\b;
%% 3.由L_转置*x=y得到方程组的解x
x=L'\y%输出线性方程组的根

function x = Cholesky_method(A,b)
%Cholesky平方根法解方程组
%A为方程组的系数矩阵 b为方程组的右端项;
    n = length(A);
    L = zeros(n);
    for k = 1:n
        delta = A(k,k);
        for j = 1:k-1
            delta = delta-L(k,j)^2;
        end

        L(k,k) = sqrt(delta);
        for i = k+1:n
            L(i,k) = A(i,k);
            for j = 1:n-1
                L(i,k) = L(i,k)-L(i,j)*L(k,j);
            end
            L(i,k) = L(i,k)/L(k,k);
        end

    end
    L 
    x =zeros(n,1);
    y =zeros(n,1);
    y(1) = b(1)/L(1,1);
    for i = 2:n
        ly = 0;
        for j = 1:i-1
            ly = ly+L(i,j)*y(i);
        end
        y(i) = (b(i)-ly)/L(i,i);
    end
    x(n) = y(n)/L(n,n);
    for i = n-1:-1:1
        lx = 0;
        for j = i+1:n
            lx = lx+L(j,i)*x(j);
        end
        x(i) = (y(i)-lx)/L(i,i);
    end
end

%Cholesky平方根法解方程组
%A为方程组的系数矩阵 b为方程组的右端项;
A6= [1 2 -1;2 5 1; -1 1 14];
b6 = [3;4;3];
x6 = Cholesky_method(A6,b6);
disp(['方程组的解：x6= ']);
disp(x6)

Cholesky分解--改进平方根法

function x=chol_ldlt_method(A,b)
%function x=chol_ldlt_method(A,b)
%Cholesky改进平方根法解方程组
%A为方程组的系数矩阵 b为方程组的右端项;
    n = length(A);
    L = eye(n);
    D = zeros(n);
    d = zeros(1,n);
    T = zeros(n);
    for k =1:n
        d(k) = A(k,k);
        for j = 1:k-1
            d(k)=d(k)-L(k,j)*T(k,j);
        end
        for i=k+1:n
            T(i,k) = A(i,k);
            for j = 1:k-1
                T(i,k) =T(i,k) -T(i,j)*L(k,j);
            end
            L(i,k) = T(i,k)/d(k);
        end
    end
    D = diag(d);
    L
    D
    x =zeros(n,1);
    y =zeros(n,1);
    d1 = zeros(n,1);
    d1 = diag(D);
    y(1) = b(1);
    for i =2:n
        ly = 0;
        for k=1:i-1
            ly = ly+L(i,k)*y(k);
        end
        y(i) = b(i)-ly;
    end
    x(n) = y(n)/d1(n);
    for i = n-1:-1:1
        lx = 0;
        for k=i+1:n
            lx = lx+L(k,i)*x(k);
        end
        x(i) = y(i)/d1(i)-lx;
end
%function x=chol_ldlt_method(A,b)
%Cholesky改进平方根法解方程组
%A为方程组的系数矩阵 b为方程组的右端项;
A7= [1 2 -1;2 5 1; -1 1 14];
b7 = [3;4;3];
x7=chol_ldlt_method(A7,b7);
disp(['方程组的解：x7= ']);
disp(x7)

奇异值分解法

%奇异值分解计算线性方程组
a=[6.5 -1 -1 3.6
6.2 7 -5 4
3 2.1 -6 4.8
1 5.6 3.7 2.1];
b=[12.3 21.4 -7.8 21]';
[u,s,v]=svd(a)
x=v*inv(s)*u'*b