前言

上一篇文章是我第一篇阅读量破百的Matlab文章，有点小鸡冻，就把它发给了当时教我相关内容的老师，收获了不小鼓励（自吹一下哈，不要脸了属于是）。然后晚上打开文章自我欣赏的时候，突然意识到，不对劲，我居然光巴巴讲了半天理论，没讲应用。这就像发图不留*，是一种被人唾弃的行为。

所以本节来讲点实用的。讲讲Matlab里面的“内置”函数。

一、上一节遗留的问题——lu函数

上一节从数学角度分析了为什么要进行LU分解，以及怎么进行LU分解。了解了原理，就可以将之运用在代码中就可以找到相应的内置函数。

在冲动编函数之前，翻了翻书，果然，Matlab给了一个内置函数“lutx”。当我兴冲冲打开编辑器准备查看时，居然提醒我不是内置函数？？？

作为一个老低血压人士，我绷住了。然后一番实验，终于找到了那个内置函数——“lu”。看来应该是我的Matlab版本比较新，书比较旧。因为版本更替而更换内置函数在Matlab里面很常见，所以不推荐用非常新版本的Matlab。我用的是2018.

function [L,U,x]=lu(A,b)
[n,n]=size(A);
p=eye(n);
for k=1:n-1[r,m]=max(abs(A(k:n,k)));  %选择主元r存值，m存行m=m+k-1;if(A(m,k)~=0)  %  确认这个主元不为0，否侧会产生奇异解if(m~=k)   %  如果m和k行不一样，就交换A([k m],:)=A([m k],:);p([k m])=p([m k]);endfor i=k+1:n    %  前向消去步骤A(i,k)=A(i,k)/A(k,k);  %  计算乘子j=k+1:n;A(i,j)=A(i,j)-A(i,k)*A(k,j);   %  把计算后的U矩阵存在A的上三角部分endend
end
L=tril(A,-1)+eye(n,n); % L矩阵是主轴以下的部分，需要补上一个单位矩阵
U=triu(A);

如上所示，lu函数先选主元，然后进行交换，防止产生舍入误差。然后再对k行主元下方的元素进行消零。然后用A矩阵下半部分保存每一行产生的乘子，用含主轴的上半部分保存每一行运算过后留下的值。最后将其分开，分别是L和U矩阵。

值得说的是，这个内置函数的代码凝炼度非常高，初学者建议自行编写，而后与内置函数进行比对。

二、完整的求解函数——bslashtx函数

还是书上说，Matlab中内置了一个bslashtx函数，是其反斜线运算符的简化版本（就是前面说的那个左除）。

然鹅！我的2018版本Matlab还是找不到这个函数！！！（是不是我下Matlab的时候丢包了？？）

然后去官网搜刮了下，找到了全文代码，补充进我的Matlab函数了，代码如下：

function x = bslashtx(A,b)
% BSLASHTX Solve linear system (backslash)
% x = bslashtx(A,b) solves A*x = b
[~,n] = size(A);
if isequal(triu(A,1),zeros(n,n))% Lower triangularx = forward(A,b);return
elseif isequal(tril(A,-1),zeros(n,n))% Upper triangularx = backsubs(A,b);return
elseif isequal(A,A')[R,fail] = chol(A);if ~fail% Positive definitey = forward(R',b);x = backsubs(R,y);returnend
end
% Triangular factorization
[L,U,p] = lu(A);
% Permutation and forward elimination
y = forward(L,p*b);
% Back substitution
x = backsubs(U,y);
% ------------------------------
function x = forward(L,x)
% FORWARD. Forward elimination.
% For lower triangular L, x = forward(L,b) solves L*x = b.
[~,n] = size(L);
x(1) = x(1)/L(1,1);
for k = 2:nj = 1:k-1;x(k) = (x(k) - L(k,j)*x(j))/L(k,k);
end
% ------------------------------
function x = backsubs(U,x)
% BACKSUBS. Back substitution.
% For upper triangular U, x = backsubs(U,b) solves U*x = b.
[~,n] = size(U);
x(n) = x(n)/U(n,n);
for k = n-1:-1:1j = k+1:n;x(k) = (x(k) - U(k,j)*x(j))/U(k,k);
end

这个“内置”函数包括了三个部分。

第一部分是判断所输入系数矩阵的形状。总共有三种特殊类型：上三角矩阵（Upper triangular）、下三角矩阵（Lower triangular）和正定矩阵（Positive definite）。如果发现系数矩阵是这三种之一，可以直接快速解开，不需要进行LU分解等变换。

第二部分是进行lu分解。直接调用lu函数，将A矩阵分解成L、U和p，这个过程参考上一节和本节的上一小节。

第三部分是该函数内部自己设置的两个函数。分别对应于求系数矩阵为上三角矩阵和下三角矩阵。这段代码在上上一节讲过。可以看到当时给出的是求上三角矩阵的两个例子，求下三角矩阵时，只需要将循环方向改变即可。而这种方向性的改变不能统一归为一个函数，所以必须把上三角和下三角分开讨论。

三、一个好玩的函数——lugui函数

这也算是一个“内置”函数——对没错是我电脑里找不到的“内置”函数，也是在官网扒了代码。相信很多和我用同一版本的人也没有这些代码，所以把他们都摘出来了。

function [L,U,pv,qv] = lugui(A,pivotstrat)
%LUGUI Gaussian elimination demonstration.
%
%   LUGUI(A) shows the steps in LU decomposition by Gaussian elimination.
%   At each step of the elimination, the pivot that would be chosen by
%   MATLAB's partial pivoting algorithm is shown in magenta.  You can use
%   the mouse to pick any pivot.  The pivot is shown in red, the emerging
%   columns of L in green, and the emerging rows of U in blue.
%
%   LUGUI with no arguments uses a random integer test matrix.
%   Type 'help golub' for a description of the test matrices.
%
%   A popup menu allows the pivot strategy to be changed dynamically.
%   lugui(A,'pick'), choose pivots with the mouse.
%   lugui(A,'diagonal'), use diagonal elements as pivots.
%   lugui(A,'partial'), use the largest element in the current column.
%   lugui(A,'complete'), use the largest element in the unreduced matrix.
%
%   [L,U,p,q] = lugui(A,...) returns a lower triangular L, an upper
%   triangular U and permutation vectors p and q so that L*U = A(p,q).
%
%   See also PIVOTGOLF.% Initializeif nargin < 2pivotstrat = 'pick';
end
if nargin < 1n = 2 + ceil(6*rand);A = golub(n);
end
Asave = A;[m,n] = size(A);
shg
clf
dx = 100;
dy = 30;
warns = warning('off','MATLAB:divideByZero');
set(gcf,'double','on','name','LU Gui', ...'menu','none','numbertitle','off','color','white', ...'pos',[480-(dx/2)*min(9,n) 320 (n+1)*dx (m+3)*dy], ...'windowbuttonupfcn','set(gcf,''tag'',''pivot'')')
stop = uicontrol('style','toggle','string','X','fontweight','bold', ...'back','w','pos',[(n+1)*dx-25 (m+3)*dy-25 25 25]);
axes('pos',[0 0 1 1])
axis off
Lcolor = [0 .65 0];
Ucolor = [0  0 .90];
Acolor = [0 0 0];
PartialPivotColor = [1 0 1];
PivotColor = [1 0 0];
TempColor = [1 1 1];
paws = 0.1;% Each element has its own handlefor j = 1:nfor i = 1:mt(i,j) = text('units','pixels','string',spf(A(i,j)), ...'fontname','courier','fontweight','bold','fontsize',14, ...'horiz','right','color',Acolor, ...'pos',[20+j*dx 20+(m+2-i)*dy],'userdata',[i j], ...'buttondownfcn','set(gcf,''userdata'',get(gco,''userdata''))');end
end% Menusswitch lower(pivotstrat)case 'pick', val = 1;case 'diagonal', val = 2;case 'partial', val = 3;case 'complete', val = 4;otherwise, val = 1;
end
pivotstrat = uicontrol('pos',[60+(dx/2)*(n-2) 20 180 20],'style','pop', ...'val',val,'fontsize',12,'back','white','string',{'Pick a pivot', ...'Diagonal pivoting','Partial pivoting','Complete pivoting'});% Eliminationpv = 1:m;
qv = 1:n;
for k = 1:min(m,n)%  If possible, quit earlyif all(all(A(k:m,k:n)==0)) || all(all(~isfinite(A(k:m,k:n))))for l = k:min(m,n)for i = l+1:mset(t(i,l),'string',spf(A(i,l)),'color',Lcolor)drawnowendfor j = l:nset(t(l,j),'string',spf(A(l,j)),'color',Ucolor)drawnowendendbreakendif (m == n) && (k == n)p = n;q = n;elsepp = min(find(abs(A(k:m,k)) == max(abs(A(k:m,k)))))+k-1;set(t(pp,k),'color',PartialPivotColor)p = 0;q = 0;while p < k || q < k || p > m || q > nswitch get(pivotstrat,'val')case 1  % Pick a pivot with mousepq = get(gcf,'userdata');if isequal(get(gcf,'tag'),'pivot') && ~isempty(pq)p = pq(1);q = pq(2);set(gcf,'tag','','userdata',[])elsedrawnowendcase 2 % Diagonal pivotingp = k;q = k;case 3 % Partial pivotingp = pp;q = k;case 4 % Complete pivoting[p,q] = find(abs(A(k:m,k:n)) == max(max(abs(A(k:m,k:n)))));p = p(1)+k-1;q = q(1)+k-1;endif get(stop,'value') == 1, break, endendif get(stop,'value') == 1, break, endset(t(pp,k),'color',Acolor)set(t(p,q),'color',PivotColor)endif get(stop,'value') == 1, break, endpause(10*paws)%  Swap columnsA(:,[q,k]) = A(:,[k,q]);qv([q,k]) = qv([k,q]);for s = .05:.05:1for i = 1:mset(t(i,k),'pos',[20+(k+s*(q-k))*dx 20+(m+2-i)*dy])set(t(i,q),'pos',[20+(q+s*(k-q))*dx 20+(m+2-i)*dy])enddrawnowendt(:,[q,k]) = t(:,[k,q]);for i = 1:mset(t(i,k),'string',spf(A(i,k)),'userdata',[i k])set(t(i,q),'string',spf(A(i,q)),'userdata',[i q])endpause(10*paws)%  Swap rowsA([p,k],:) = A([k,p],:);pv([p,k]) = pv([k,p]);for s = .05:.05:1for j = 1:nset(t(k,j),'pos',[20+j*dx 20+(m+2-(k+s*(p-k)))*dy])set(t(p,j),'pos',[20+j*dx 20+(m+2-(p+s*(k-p)))*dy])enddrawnowendt([p,k],:) = t([k,p],:);pause(10*paws)for j = k:nset(t(k,j),'string',spf(A(k,j)),'userdata',[k j])set(t(p,j),'string',spf(A(p,j)),'userdata',[p j])endpause(10*paws)%  Skip step if column is all zeroif all(A(k:m,k) == 0)for i = k+1:mset(t(i,k),'string',spf(A(i,k)),'color',Lcolor)drawnowendfor j = k:nset(t(k,j),'string',spf(A(k,j)),'color',Ucolor)drawnowendelse%     Compute multipliers in Lfor i = k+1:mA(i,k) = A(i,k)/A(k,k);set(t(i,k),'string',spf(A(i,k)),'color',Lcolor)pause(paws)drawnowend%     Elimination    for j = k+1:nfor i = k+1:mset(t(i,j),'color',TempColor)drawnowpause(paws)A(i,j) = A(i,j) - A(i,k)*A(k,j);set(t(i,j),'string',spf(A(i,j)),'color',Acolor)drawnowpause(paws)endendfor j = k:nset(t(k,j),'string',spf(A(k,j)),'color',Ucolor)drawnowendpause(paws)endif k < min(m,n), pause(10*paws), end
end% Seperate L and U into two matricesdelete(pivotstrat)for s = .1:.1:1.5for j = 1:nfor i = 1:mif i <= jset(t(i,j),'pos',[20+(j+.10*s)*dx 20+(m+2-i)*dy])elseset(t(i,j),'pos',[20+(j-.10*s)*dx 20+(m+2-s-i)*dy])endendenddrawnow
end% Insert ones on diagonal of Lr = min(m,n);
for j = 1:rtext('units','pixels','string',spf(1.0), ...'fontname','courier','fontweight','bold','fontsize',14, ...'horiz','right','color',Lcolor, ...'pos',[20+(j-0.15)*dx 20+(m+.5-j)*dy]);
end
drawnow
warning(warns)if nargout > 0L = tril(A(:,1:r),-1) + eye(m,r);U = triu(A(1:r,:));
elseset(gcf,'userdata',Asave)set(stop,'value',0,'callback','close(gcf)')uicontrol('pos',[(n+1)*dx-70 10 60 20],'string','repeat', ...'back','w','fontsize',12,'callback','lugui(get(gcf,''userdata''))')
end%------------------------------------------------------------function s = spf(aij)
% Subfunction to format text strings
if aij == 0f = '%10.0f';
elseif (abs(aij) < 1.e-4) || (abs(aij) >= 1.e4) f = '%10.1e';
elsef = '%10.4f';
end
s = sprintf(f,aij);

这个代码颇具复杂性，初学者不需要掌握，看懂其中某些步骤即可。它使用了Matlab的GUI功能，可以显示高斯消元法的每个步骤。可以用于讲课和考研题目的演算（bushi）。

比如我还是输入之前的那个系数矩阵A=[10 -7 0;-3 2 6;5 -1 5]; 调用该函数lugui(A)

然后出现如下图的界面：

这图中，粉红色的元素是系统推荐的主元。在这个例子中，我们在处理第二行的时候，进行了行变换，因为主元在第三行。可以通过这个gui来验证一下。首先点击10.0000，在手动模式下，点击选中的主元会变成红色。

而后，程序完成第一列的消元处理，第一行被固定，变成蓝色。第一列也被固定，变成绿色。由上一小节的内容可以知道，下方存储的是乘子矩阵L，上方存储的是消元后的矩阵U。同时，返回了第二个系统推荐的主元。

再次选择第三行第二列的元素作为我们想要的主元。然后它变成红色。同时完成行变换以及消元。

最后，通过判断行数，得出矩阵已经LU分解完毕，将矩阵拆开，返回上方的U矩阵和下方的L矩阵。按repeat可以返回重新玩一遍。

对于主元选择功能可以有四种选择：

1、手工选主元：就是上面讲的那样，用鼠标来点击想要的主元；

2、对角线选主元：默认用对角线元素作为主元；

3、部分选主元：自动使用lu函数的策略来选主元；

4、完全选主元：用剩余子阵里绝对值最大的元素作为主元。

总结

至此，已经学会了用Matlab去求解线性方程。所用的方法是高斯消元法，使用了LU分解这种矩阵思想。这对于一些工程应用基础和线性代数答案检查都有帮助。

内置的lu、bslashtx、lugui函数——Matlab解线性方程组（4）相关推荐

未捕获的错误：始终违反：元素类型无效：预期为字符串（对于内置组件）或类/函数，但得到了：对象
本文翻译自:Uncaught Error: Invariant Violation: Element type is invalid: expected a string (for built-in ...
python turtle库setpos_Python内置海龟(turtle)库绘图命令详解（二）
继续谈利用海龟库(turtle库)做图.在这篇文章(Python内置海龟(turtle)库绘图命令详解(一))中已经介绍了turtle的一些基本画图命令,包括画布的设计.画笔属性与状态的设置以及画笔的 ...
matlab在解线性方程组的应用,matlab解线性方程组线性方程组及MATLAB应用
matlab解线性方程组线性方程组及MATLAB应用 1matlab 解线性方程组线性方程组及 MATLAB 应用数值实验线性方程组与 MATLAB 应用王1．实验目的:理解矩阵的范数与条件数. ...
laravel框架内置的各种路径帮助函数
一.应用场景: 当需要引入第三方路径,而通过传统的'./'或者'../'无法准确定位的时候,我们需要能直接获取项目所在网站的根目录,然后根据根目录找到我们需要引入的第三方包. 二.例子首先声明,博主 ...
mysql5720_Mysql内置功能《五》函数
一函数 MySQL中提供了许多内置函数,例如: 一.数学函数 ROUND(x,y) 返回参数x的四舍五入的有y位小数的值 RAND() 返回0到1内的随机值,可以通过提供一个参数(种子)使RAND( ...
python序列类型-Python内置序列类型之集合类型详解
1.集合概念具有某种特定性质的事物的总体,集合里的东西叫作元素.Python中,集合(set)是一个无序不重复元素的序列. 2.集合的创建可以使用大括号 { } 或者 set() 函数创建集合,注 ...
Dubbo内置4种负载均衡算法(详解)
1.1 什么是负载均衡在实际开发中,一个服务基本都是集群模式的,也就是多个功能相同的项目在运行,这样才能承受更高的并发,这时一个请求到这个服务,就需要确定访问哪一个服务器 Dubbo框架内部支持负载 ...
matlab解线性方程组后结果是小数,MATLAB线性方程组求解
有唯一解线性方程组求法对于一般的,有唯一解的线性方程组,我们可以转换成矩阵的形式: A x = b Ax=bAx=b 则可以用矩阵运算求解x,即x=A\b 有无穷解的线性方程组求法齐次线性方程组的 ...
MATLAB解线性方程组
rref 函数把矩阵换为行最简形可以用来解线性方程组,求矩阵的秩,求矩阵行最简形(每行首元所在的列只有它一个是1)首元所在的列数. 例如我们知道一个方程组 A*X=b 中 A 系数矩阵和b ...

内置的lu、bslashtx、lugui函数——Matlab解线性方程组（4）

前言

一、上一节遗留的问题——lu函数

二、完整的求解函数——bslashtx函数

三、一个好玩的函数——lugui函数

总结

内置的lu、bslashtx、lugui函数——Matlab解线性方程组（4）相关推荐

最新文章

热门文章