线性方程组求解

胡枝子.yue

已于 2022-10-22 10:03:03 修改

阅读量470

点赞数

分类专栏：数值计算文章标签： matlab

于 2022-10-22 10:01:50 首次发布

本文链接：https://blog.csdn.net/weixin_52508924/article/details/127457717

版权

数值计算专栏收录该内容

3 篇文章 0 订阅

订阅专栏

数值计算实验的简单记录

文章目录

一、列主元消元法
二、Guass-Jordan消元法
三、Guass-Jordan法求逆矩阵
四、平方根法
五、改进平方根法
六、LU分解

前言

各算法的核心思想可查阅金一庆的《数值方法》第二版。如果高代学得不错，也可以直接百度算法原理来看，丘维声高代中基本都有提过。

一、列主元消元法

A=input('请输入系数矩阵：')
%A=[2 -1 3;-1 -2 3;1 3 1]
b=input('请输入非齐次向量b：')
%b=[4 5 6]
b=b'
B=[A,b]
[h,l]=size(A)
[h1,l1]=size(b)
if h~=h1
    disp('A或b存在输入错误！')
elseif rank(A)~=rank(B)
    disp('该线性方程无解！')
elseif det(A)~=0
    disp('该方程存在唯一解。')
end
%%若方程存在唯一解，我们就可以进一步使用列主元消元法求解方程
for i=1:h
    ma=max(abs(B(i:end,i)))
    [x,y]=find(B==ma)%[x,y]=find(max(abs(B(1:end,1))))
    B([i,x],:)=B([x,i],:)%B([1,x],:)=B([x,1],:)
    for j=i+1:h
        B(j,:)=B(j,:)-B(i,:).*B(j,i)/B(i,i)%B(2,:)=B(2,:)-B(1,:).*(B(2,1)/B(1,1))
    end
end
for k=1:h
    for j1=1:h-k
        B(j1,:)=B(j1,:)-B(h+1-k,:)*B(j1,h+1-k)/B(h+1-k,h+1-k)
%B(2,:)=B(2,:)-B(h+1-1,:)*B(2,h+1-1)/B(h+1-1,h+1-1)
    end
end
xx=B(:,end)./diag(B(:,1:h));
xx=xx';
t=['该线性方程组的解为：',num2str(xx)];
disp(t)
%inv(A)*b

二、Guass-Jordan消元法

%A=input('请输入系数矩阵：')
A=[2 -3 1;-3 5 -1;1 -1 2]
%b=input('请输入非齐次向量b：')
b=[3 -6 -2]
b=b'
B=[A,b]
[h,l]=size(A)
[h1,l1]=size(b)
if h~=h1
    disp('A或b存在输入错误！')
elseif rank(A)~=rank(B)
    disp('该线性方程无解！')
elseif det(A)~=0
    disp('该方程存在唯一解。')
end
%%若方程存在唯一解，我们就可以进一步使用列主元消元法求解方程
for i=1:h
    if B(i,i)==0
        [x,y]=find(A(:,i)~=0)
        B([i,x],:)=B([x,i],:)
    end
    B(i,:)=B(i,:)/B(i,i)
     for j=i+1:h
        B(j,:)=B(j,:)-B(i,:).*B(j,i)/B(i,i)%B(2,:)=B(2,:)-B(1,:).*(B(2,1)/B(1,1))
    end
end
for k=1:h
    for j1=1:h-k
        B(j1,:)=B(j1,:)-B(h+1-k,:)*B(j1,h+1-k)/B(h+1-k,h+1-k)
%B(2,:)=B(2,:)-B(h+1-1,:)*B(2,h+1-1)/B(h+1-1,h+1-1)
    end
end
t=['该线性方程组的解为：',num2str((B(:,end))')];
disp(t)

三、 Guass-Jordan法求逆矩阵

%A=[-3 8 5;2 -7 4;1 9 -6]
%A=[2 1 -3 -1;3 1 0 7;-1 2 4 -2;1 0 -1 5]
[h,l]=size(A)
B=[A,eye(h)]
if h==l && det(A)
    disp('矩阵A可逆')
else
    disp('矩阵没有逆矩阵')
end
for i=1:h
    if B(i,i)==0
        [x,y]=find(A(:,i)~=0)
        B([i,x],:)=B([x,i],:)
    end
    B(i,:)=B(i,:)/B(i,i)
     for j=i+1:h
        B(j,:)=B(j,:)-B(i,:).*B(j,i)/B(i,i)%B(2,:)=B(2,:)-B(1,:).*(B(2,1)/B(1,1))
    end
end
for k=1:h
    for j1=1:h-k
        B(j1,:)=B(j1,:)-B(h+1-k,:)*B(j1,h+1-k)/B(h+1-k,h+1-k)
%B(2,:)=B(2,:)-B(h+1-1,:)*B(2,h+1-1)/B(h+1-1,h+1-1)
    end
end
A_inv=B(:,h+1:end)

四、平方根法

% A=input('请输入系数矩阵A：')
tic
t1=clock
%A=[6,1,0;1,4,1;0,1,14]
%b=input('请输入非齐次向量b：')
%b=[6 24 322]
 A=[2 -1 3;-1 -2 3;1 3 1]
 b=[4 5 6]
b=b'
[h,l]=size(A)
if A~=A'  
    disp('该方程不是对称矩阵')
elseif size(A,1)==sum(eig(A)>0)
    disp('该方程是对称正定阵，可以使用cholesky分解！')
else
    disp('该方程是对称矩阵，但不正定')
end
%判断矩阵是否可以cholesky分解
L=zeros(h)
for i=1:h
    L(i,i)=sqrt(A(i,i)-sum(L(i,1:i-1).*L(i,1:i-1)))
    for j=i+1:h
    L(j,i)=(A(j,i)-sum(L(j,1:i-1).*L(i,1:i-1)))/L(i,i)
    end
end
C=L
L=[C,b]
for k=1:h
    for j1=k+1:h
        L(j1,:)=L(j1,:)-L(k,:)*L(j1,k)/L(k,k)
    end
end
y=L(:,end)./diag(L(:,1:h));
D=[C',y]
for k1=1:h
    for j1=1:h-k1
        D(j1,:)=D(j1,:)-D(h+1-k1,:)*D(j1,h+1-k1)/D(h+1-k1,h+1-k1)
%B(2,:)=B(2,:)-B(h+1-1,:)*B(2,h+1-1)/B(h+1-1,h+1-1)
    end
end
toc
disp(['toc计算最后一次循环运行时间',num2str(toc)])
disp(['etime程序总运行时间：',num2str(etime(clock,t1))])
L=C
x=D(:,end)./diag(D(:,1:h))

五、改进平方根法

% A=input('请输入系数矩阵A：')
tic
t1=clock;
% A=[6,1,0;1,4,1;0,1,14]
% b=[6 24 322]
%b=input('请输入非齐次向量b：')
 A=[2 -1 3;-1 -2 3;1 3 1]
 b=[4 5 6]
b=b'
[h,l]=size(A)
if A~=A'  
    disp('该方程不是对称矩阵')
elseif size(A,1)==sum(eig(A)>0)
    disp('该方程是对称正定阵，可以使用cholesky分解！')
else
    disp('该方程是对称矩阵，但不正定')
end
%判断矩阵是否可以cholesky分解
d=zeros(h)
d(1,1)=A(1,1)
T=zeros(h)
for i=2:h
    T(i,1)=A(i,1)
    l(i,1)=T(i,1)/d(1,1)
    l(1,1)=1
    %d(2,2)=A(2,2)-sum(T(2,1).*l(2,1))
end
for i=2:h
    d(i,i)=A(i,i)-sum(T(i,i-1).*l(i,i-1))
    for j=2:i
    T(i,j)=A(i,j)-sum(T(i,1:j-1).*l(j,1:j-1))
    T(i,i)=0
    l(i,j)=T(i,j)/d(j,j)
    l(i,i)=1
    
    end     
end

L=[l*d,b]
for k=1:h
    for j1=k+1:h
        L(j1,:)=L(j1,:)-L(k,:)*L(j1,k)/L(k,k)
    end
end
y=L(:,end)./diag(L(:,1:h));
D=[l',y]
for k1=1:h
    for j1=1:h-k1
        D(j1,:)=D(j1,:)-D(h+1-k1,:)*D(j1,h+1-k1)/D(h+1-k1,h+1-k1)
%B(2,:)=B(2,:)-B(h+1-1,:)*B(2,h+1-1)/B(h+1-1,h+1-1)
    end
end
toc
disp(['toc计算最后一次循环运行时间',num2str(toc)])
disp(['etime程序总运行时间：',num2str(etime(clock,t1))])

x=D(:,end)./diag(D(:,1:h))

六、LU分解

% A=input('请输入系数矩阵A：')
A=[2 -1 3;-1 -2 3;1 3 1]
[h,l]=size(A)
for i=1:h-1
    Lt=inv(L)
    L=eye(h)
    for j=i+1:h
        L(j,i)=-A(j,i)/A(i,i)
    end
    Lt=Lt*inv(L)
    A=L*A
end
U=A
L=Lt
%lu(A)