目录
前言
上一篇文章是我第一篇阅读量破百的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]);
end
for 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的上三角部分
end
end
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 triangular
x = forward(A,b);
return
elseif isequal(tril(A,-1),zeros(n,n))
% Upper triangular
x = backsubs(A,b);
return
elseif isequal(A,A')
[R,fail] = chol(A);
if ~fail
% Positive definite
y = forward(R',b);
x = backsubs(R,y);
return
end
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:n
j = 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:1
j = 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.
% Initialize
if nargin < 2
pivotstrat = 'pick';
end
if nargin < 1
n = 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 handle
for j = 1:n
for i = 1:m
t(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
% Menus
switch 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'});
% Elimination
pv = 1:m;
qv = 1:n;
for k = 1:min(m,n)
% If possible, quit early
if 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:m
set(t(i,l),'string',spf(A(i,l)),'color',Lcolor)
drawnow
end
for j = l:n
set(t(l,j),'string',spf(A(l,j)),'color',Ucolor)
drawnow
end
end
break
end
if (m == n) && (k == n)
p = n;
q = n;
else
pp = 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 > n
switch get(pivotstrat,'val')
case 1 % Pick a pivot with mouse
pq = get(gcf,'userdata');
if isequal(get(gcf,'tag'),'pivot') && ~isempty(pq)
p = pq(1);
q = pq(2);
set(gcf,'tag','','userdata',[])
else
drawnow
end
case 2 % Diagonal pivoting
p = k;
q = k;
case 3 % Partial pivoting
p = 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;
end
if get(stop,'value') == 1, break, end
end
if get(stop,'value') == 1, break, end
set(t(pp,k),'color',Acolor)
set(t(p,q),'color',PivotColor)
end
if get(stop,'value') == 1, break, end
pause(10*paws)
% Swap columns
A(:,[q,k]) = A(:,[k,q]);
qv([q,k]) = qv([k,q]);
for s = .05:.05:1
for i = 1:m
set(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])
end
drawnow
end
t(:,[q,k]) = t(:,[k,q]);
for i = 1:m
set(t(i,k),'string',spf(A(i,k)),'userdata',[i k])
set(t(i,q),'string',spf(A(i,q)),'userdata',[i q])
end
pause(10*paws)
% Swap rows
A([p,k],:) = A([k,p],:);
pv([p,k]) = pv([k,p]);
for s = .05:.05:1
for j = 1:n
set(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])
end
drawnow
end
t([p,k],:) = t([k,p],:);
pause(10*paws)
for j = k:n
set(t(k,j),'string',spf(A(k,j)),'userdata',[k j])
set(t(p,j),'string',spf(A(p,j)),'userdata',[p j])
end
pause(10*paws)
% Skip step if column is all zero
if all(A(k:m,k) == 0)
for i = k+1:m
set(t(i,k),'string',spf(A(i,k)),'color',Lcolor)
drawnow
end
for j = k:n
set(t(k,j),'string',spf(A(k,j)),'color',Ucolor)
drawnow
end
else
% Compute multipliers in L
for i = k+1:m
A(i,k) = A(i,k)/A(k,k);
set(t(i,k),'string',spf(A(i,k)),'color',Lcolor)
pause(paws)
drawnow
end
% Elimination
for j = k+1:n
for i = k+1:m
set(t(i,j),'color',TempColor)
drawnow
pause(paws)
A(i,j) = A(i,j) - A(i,k)*A(k,j);
set(t(i,j),'string',spf(A(i,j)),'color',Acolor)
drawnow
pause(paws)
end
end
for j = k:n
set(t(k,j),'string',spf(A(k,j)),'color',Ucolor)
drawnow
end
pause(paws)
end
if k < min(m,n), pause(10*paws), end
end
% Seperate L and U into two matrices
delete(pivotstrat)
for s = .1:.1:1.5
for j = 1:n
for i = 1:m
if i <= j
set(t(i,j),'pos',[20+(j+.10*s)*dx 20+(m+2-i)*dy])
else
set(t(i,j),'pos',[20+(j-.10*s)*dx 20+(m+2-s-i)*dy])
end
end
end
drawnow
end
% Insert ones on diagonal of L
r = min(m,n);
for j = 1:r
text('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 > 0
L = tril(A(:,1:r),-1) + eye(m,r);
U = triu(A(1:r,:));
else
set(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 == 0
f = '%10.0f';
elseif (abs(aij) < 1.e-4) || (abs(aij) >= 1.e4)
f = '%10.1e';
else
f = '%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分解这种矩阵思想。这对于一些工程应用基础和线性代数答案检查都有帮助。
7535
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
