-
定解方程组
计算过程:
注意:
矩阵系数中,逗号和分号的使用
矩阵除法用反斜杠\
矩阵中的值提取用小括号,行列数从1开始算
-
不定方程组
未知数大于方程数,其解有无数个,matlab可以求出其中一个特定解
计算过程
以三元二次方程组为例(其他情况可类推):
{
17
x
2
+
9
y
2
+
5
z
2
−
24
x
y
+
12
y
z
−
14
x
z
−
6
x
−
6
z
=
0
x
+
2
y
−
z
+
1
=
0
\begin{cases} 17x^2+9y^2+5z^2-24xy+12yz-14xz-6x-6z=0\\ x+2y-z+1=0 \end{cases}
{17x2+9y2+5z2−24xy+12yz−14xz−6x−6z=0x+2y−z+1=0
在Matlab命令行窗口输入以下命令:
syms x y z
eq1=17*x^2+9*y^2+5*z^2-24*x*y+12*y*z-14*x*z-6*x-6*z
eq2=x+2*y-z+1
[x,y,z]=solve(eq1,eq2,x,y,z)
得到以下结果:
x =
(19^(1/2)*6i)/125 - 21/125
- (19^(1/2)*6i)/125 - 21/125
y =
- (19^(1/2)*3i)/125 - 52/125
(19^(1/2)*3i)/125 - 52/125
z =
0
0
计算过程:
注意:
矩阵系数中,逗号和分号的使用
矩阵除法用反斜杠\
矩阵中的值提取用小括号,行列数从1开始算
未知数大于方程数,其解有无数个,matlab可以求出其中一个特定解
计算过程
转载于:https://www.cnblogs.com/derek32/p/4042829.html
综合实例应用:方程组的求解
无论工程应用问题,还是数学计算问题,方程组都是解决问题转化的重要途径之一,将复杂问题转化为简单的方程组矩阵求解问题。
>> %创建方程组系数矩阵
>> 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
利用矩阵分解来求解线性方程组,是工程计算中最常用的计算。
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
利用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
综合实例—方程组的求解,到这里就结束啦!感谢观看,希望这篇文章对大家有帮助。
前面有两篇博文分别介绍了:
MATLAB求常微分方程的解析解
MATLAB求常微分方程的数值解
为了形成一个体系,我决定把普通方程组的求解也介绍一下。
本博文也是按照MATLAB的官方文档展开的(推荐大家多看官方文档)
一般形式 S=solve(eqns,vars,Name,Value)
,其中:
eqns是需要求解的方程组;
vars是需要求解的变量;
Name-Value对用于指定求解的属性(一般用不到);
S是结果,对应于vars中变量;
方程:sin(x)=1
代码:
syms x; %定义x是一个未知量
eqn=sin(x)==1; % 定义方程,eqn只是一个代号,代表sin(x)==1
solX=solve(eqn,x) % 求方程eqn中的x,放入solX中
结果:
说明: MATLAB定义方程用的是 == 符号,就是这样规定的哈。
注意: 细心的同学应该发现了,本例的解实际上应该是 pi/2+2k*pi ,怎么得到呢?
添加Name-Value对即可解决,输入以下代码:
syms x; %定义x是一个未知量
eqn=sin(x)==1; % 定义方程,eqn只是一个代号,代表sin(x)==1
[solX,params,cond]=solve(eqn,x,'ReturnConditions',true) % 求方程eqn中x的所有解,放入solX中,params是参数,cond存储参数性质
得到理想结果:
方程: ax²+bx+c=0
代码:
syms x a b c; %定义x a b c是未知量
eqn=a*x^2+b*x+c==0;% 定义方程
solX=solve(eqn,x) % 解方程
结果:
说明: 这里就简单的把未知参数用syms声明就可以了。
方程:
代码:
syms u v; % 定义u v 是未知量
eqns=[2*u+v==0,u-v==1]; % 定义方程组
vars=[u,v]; % 定义求解的未知量
[solU,solV]=solve(eqns,vars) % 求解eqns中的vars未知量,分别存储
sol=solve(eqns,vars); % 求解eqns中的vars未知量,以结构体的形式存储到sol中
solU1=sol.u % 从sol结构体中取出变量u的解
solV1=sol.v % 从sol结构体中取出变量v的解
结果:
说明: 本例中有两个求解的变量,有两种存储方式,已在代码中介绍。
方程: sin(x)==x²-1
代码:
syms x; % 定义x是未知量
fplot(sin(x),[-2,2]); % 绘制y=sin(x)的图像
hold on;
fplot(x^2-1,[-2,2]); % 绘制y=x^2-1的图像
hold off;
eqn=sin(x)==x^2-1; % 定义方程
solX=solve(eqn,x) % 直接求解,返回其找到的第一个数值近似解
solX1=vpasolve(eqn,x,[0,2]) % vpa求解,返回其在范围[0,2]内找到的第一个数值近似解
结果:
说明: 此例中无法求得精确解,slove会返回求得的第一个数值近似解,vpasolve可以返回指定范围内第一个近似解
方程:
代码:
syms x; % 定义x是未知量
eqn=[3*x+2==0,3*x+1==0]; % 定义函数
solX=solve(eqn,x) % 求解
结果: