MATLAB:方程组的求解

MATLAB:方程组的求解综合实例应用:方程组的求解无论工程应用问题,还是数学计算问题,方程组都是解决问题转化的重要途径之一,将复杂问题转化为简单的方程组矩阵求解问题。二、利用矩阵分解求解利用矩阵分解来求解线性方程组,是工程计算中最常用的计算。编写利用LU分解法求解线性方程组Ax=b的自定义函数M文件,操作方法:2.QR分解法利用QR分解法先将系数矩阵A进行QR分解A=QR,然后解Qy=b,最后解Rx=y得到原方程组的解总结综合实例—方程组的求解,到这里就结束啦!感谢观看,希望这篇文章对大家有帮助。..

系列文章目录

综合实例应用:方程组的求解



前言

无论工程应用问题,还是数学计算问题,方程组都是解决问题转化的重要途径之一,将复杂问题转化为简单的方程组矩阵求解问题。


一、求解四元一次线性方程组

>> %创建方程组系数矩阵
>> A=[2 1 -5 1;1 -3 0 -6;0 2 -1 2;1 4 -7 6];
>> b=[8 9 -5 0]';
>> %判断方程是否有解
>> %求方程组的秩
>> r=rank(4)

r =

     1

>> B=[A,b];%创建增广矩阵
>> s=rank(B)

s =

     4

>> %r=s=n(未知数)=4,则该齐次线性方程组有唯一解。
>> %利用矩阵的逆
>> x0=pinv(A)*b

x0 =

    3.0000
   -4.0000
   -1.0000
    1.0000

二、利用矩阵分解求解

利用矩阵分解来求解线性方程组,是工程计算中最常用的计算。

1.LU分解法

LU分解法是先将系数矩阵A进行LU分解,得到LU=PA,然后解Ly=Pb,最后再解Ux=y得到原方程组的解。

编写利用LU分解法求解线性方程组Ax=b的自定义函数M文件,操作方法:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

函数solvebyLU的程序,如下所示

function x=solvebyLU(A,b)
%该函数利用LU分解法求线性方程组Ax=b的解
flag=isexist(A,b); 					%调用自定义函数isexist()判断方程组解的情况
if flag==0
    disp('该方程组无解!');
    x=[];
    return;
else
    r=rank(A);
    [m,n]=size(A);
    [L,U,P]=lu(A);
    b=P*b;
        %解Ly=b
    y(1)=b(1);
    if m>1
        for i=2:m
            y(i)=b(i)-L(i,1:i-1)*y(1:i-1)';
        end
    end
    y=y';
        %解Ux=y得原方程组得一个特解
    x0(r)=y(r)/U(r,r);
    if r>1
        for i=r-1:-1:1
            x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i);
        end
    end
    x0=x0';
     if flag==1  					%若方程组有唯一解
        x=x0;
        return;
    else        						%若方程组有无穷多解
        format rat;
        Z=null(A,'r'); 				%求出对应齐次方程组的基础解系
        [mZ,nZ]=size(Z);
        x0(r+1:n)=0;
        for i=1:nZ
            t=sym(char([107 48+i]));
            k(i)=t;      				%取k=[k1,k2...,];
        end
        x=x0;         
        for i=1:nZ          
            x=x+k(i)*Z(:,i); 			%将方程组的通解表示为特解加对应齐次通解形式
        end        
    end
end

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

函数isexist()的程序,如下所示

function y=isexist(A,b)
%该函数用来判断线性方程组Ax=b的解的存在性
%若方程组无解则返回0,若有唯一解则返回1,若有无穷多解则返回Inf。
 [m,n]=size(A);
[mb,nb]=size(b);
if m~=mb
    error('输入有误!');
    return;
end
r=rank(A);
s=rank([A,b]);
if r==s &&r==n
    y=1;
elseif r==s&&r<n
    y=Inf;
else
    y=0;
end

在这里插入图片描述

命令行代码,如下所示

>> A=[2 1 -5 1;1 -3 0 -6;0 2 -1 2;1 4 -7 6];
>> b=[8 9 -5 0]';
>> x2=solvebyLU(A,b)

x2 =

       3       
      -4       
      -1       
       1       

2.QR分解法

利用QR分解法先将系数矩阵A进行QR分解A=QR,然后解Qy=b,最后解Rx=y得到原方程组的解
在这里插入图片描述

1.编写求解线性方程组Ax=b的函数solvebyQR,代码如下:

function x=solvebyQR(A,b)
%该函数利用QR分解法求线性方程组Ax=b的解
flag=isexist(A,b); 					%调用自定义函数isexist()
if flag==0
    disp('方程组无解');
    x=[];
    return;
else
    r=rank(A);
    [m,n]=size(A);
    [Q,R]=qr(A);
    b=Q'*b;
    %解Rx=b得原方程组得一个特解
    x0(r)=b(r)/R(r,r);
    if r>1
        for i=r-1:-1:1
            x0(i)=(b(i)-R(i,i+1:r)*x0(i+1:r)')/R(i,i);
        end
    end
    x0=x0';
       if flag==1  					%若方程组有唯一解
        x=x0;
        return;
    else        						%若方程组有无穷多解
        format rat;
        Z=null(A,'r'); 				%求出对应齐次方程组得基础解系
        [mZ,nZ]=size(Z);
        x0(r+1:n)=0;
        for i=1:nZ
            t=sym(char([107 48+i]));
            k(i)=t;      					%取k=[k1,...,kr];
        end
        x=x0;         
        for i=1:nZ          
            x=x+k(i)*Z(:,i); 				%将方程组的通解表示为特解加对应齐次通解形式
        end        
    end
end

总结

综合实例—方程组的求解,到这里就结束啦!感谢观看,希望这篇文章对大家有帮助。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:https://bianchenghao.cn/39308.html

(0)
编程小号编程小号

相关推荐

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注