系列文章目录
综合实例应用:方程组的求解
前言
无论工程应用问题,还是数学计算问题,方程组都是解决问题转化的重要途径之一,将复杂问题转化为简单的方程组矩阵求解问题。
一、求解四元一次线性方程组
>> %创建方程组系数矩阵
>> 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