一直很喜欢两个老爷爷,一个是MIT主讲线性代数的Gilbert Strang,另一个就是Matlab的首席数学科学家Cleve Moler。记得作者在几年前,下载了他个人主页上的两本书。一本叫《Numerical Computing with MATLAB》,另外一本叫《Experiments with MATLAB》。当时研究了一下,就没了下文,近来我又把他的《Experiments with MATLAB》翻出来看了一下,虽然网上也有中文版,但我还是勉强的啃着英文看。

        这两天看到了书中的一篇文章,讲的是MATLAB解方程组Ax=b时最常用的反斜杠“\”。现在我打算把他的这篇文章翻译一下,分享给大家,中间可能会补充一些自己的东西。

         两个数的平均值等于3,请问这两个数分别是多少?请记住你脑子里蹦出的两个数,我们会在这一章的结尾讨论这个问题。

        我们先从一个非常简单的线性方程说起,一个未知数,一个方程。

求:

很明显答案是:

 现在我们常数换成代数,试求解:

 得到:

方程有唯一解

那么如果a=0呢?这就要看b了。如果b不等于0,那么也就没有任何x能满足:

 因此,方程无解

那如果换成b=0呢,如果b=0,那么任何x都可以满足方程:

 也就是说,方程有无穷多个解。数学家们已经对“如何判断方程是否有解?”和“如果有解,那么解是否唯一?”这两个问题讨论了几个世纪了。我们马上就会看到,那些令数学家们讨论了几个世纪的问题,对于今天的工程计算依然依然继续着。

买水果问题:

        爱丽丝买了3个苹果,12个香蕉和1个香瓜,总共花了2.36美元。鲍勃买了12个苹果和2个香瓜,总共花了5.26美元。卡罗买了2个香蕉和3个香瓜,总共花了2.77美元。请问,这三种水果的单价是多少?

我们用未知数x_{1},x_{2}x_{3}分别表示苹果,香蕉和香瓜的单价。并把上述的文字描述用方程的方式分别表示出来,得到一个三元一次方程组:

 然后,我们再用一个矩阵右乘一个未知数列向量的方式来表示上述的三元一次方程组:

 这就是我们最常见的,也是最熟悉:

 其中,A代表了3x3矩阵,x和b分别是两个3x1的列向量。b也叫right hand side。

一般而言,如果你是一个对于矩阵运算比较了解的人,那么你会马上想到用A的逆矩阵A^{-1}来解方程,即:

但是,这种解法虽然能用,但并不是一个好的想法,事实上我们不需要求出A的逆A^{-1},也能求出方程的解x。

 比如说,假设你是一个对矩阵计算完全不了解的人,基于前面求解一元一次方程ax=b的解法,你可能会想到这么算:

 从线性代数本身的理论出发,很遗憾,你不能这么算,你不能除一个矩阵。但在Matlab中,你可以这也用。

我们先输入3x3的系数矩阵A,Matlab命令如下:

A = [3 12 1; 12 0 2; 0 2 3]

Matlab会输出如下结果:

接下来我们把方程组的右边b也输进去:

 b = [2.36 5.26 2.77]'

输出:

 现在,在Matlab中求解Ax=b,就可以使用反斜杠“\”操作符了:

 这就好像是让等式Ax=b左右两边都除以A,得到:

 苹果0.29美元一个,香蕉0.05美元一个,香瓜0.89美元一个。

斜杠(除号)与反斜杠(反除号):

        对于我们最开始的一次方程7x=21而言,我们可以用反斜杠的方式表示解:

x=7\setminus 21

 其实,我们也可以用斜杠写成:

x=21/7

当然这也是我们最常用的表达方式。当然,不论是正除还是反除,他们的结果一样,都是3。

对于矩阵的运算也一样:大多数情况下,我们求解的线性方程组都可以表示成

的形式。

但在某些特殊情况下,我们的线性方程组会表示成

 对于这种情况的线性方程组,我们在求解的时候就不再使用反斜杠“\”了,而是直接使用斜杠“/”

 有时候,A\bb/A被分别称为矩阵的左除和右除,不论哪一种除法,矩阵A都是分母。

奇异矩阵singular matrix:

现在我们把前面的买水果问题做一点小小的调整。假设,卡罗同样用2.77美元,买了6个苹果和1个香瓜。则,原始的矩阵A变成:

 而原来的矩阵右边b不变:

乍一眼看上去,这组方程和原来的方程比变化不大,只是A的最后一行方程发生了改变:

 如果你继续用

 来求解这个方程的话,你会发现MATLAB会做出如下提示:

 其中,-InfInf分别表示正负无穷大,这是由于在计算的过程中遇到了除0的情况。而NaN指的是“Not-a-Number”,这是因为在计算的时候遇到了无穷大的数。

        究其原因是因为,卡罗购买的水果和所花的钱和爱丽丝,鲍勃购买的水果,以及他们所花的钱,出现了冲突。

        卡罗购买的水果品种和鲍勃一样,数量只有鲍勃的一半,但他花的钱却不等于鲍勃的二分之一。 鲍勃的一半是2.63而不是现在的2.77。矩阵A种的第三行等于第二行的一半,而方程右端b的第三个数不是第二个数的一半。这种不合理性,直接导致了Ax=b的解不存在。

如果是一半呢?

        现在我们把鲍勃的花销改成跟卡罗的花销一致,也就是维持矩阵A不变,把b改成:

这样一来我们最后的两个方程就变成下面这个样子了:

而这两个方程其实可以用一个方程来表示,多出来的一个并没有提供更多的性息。

现在我们试着用反斜杠来解一下这个方程组: 

 系统报错了,说A是奇异矩阵(并不是所有奇异矩阵MATLAB都不能求解),MATLAB解不了这个方程,我们只能自己手动求解了!

在解这种方程组时,我们应该牢记,对于奇异矩阵A而言,他的全部解等于一个特解加上一个任意实数乘以零空间

 上图中的笔记有误,把“最终得到:Ax=b的解”改成“最终得到:Ax=b的解”,同时,在上面的求解步骤中,我还使用了Matlab的另外一个函数null()来求解A的零空间也就是Ax=0的解

Matlab code:

close all
clear all
% J27,2022/08/22
%% BackSlash
%Alice buys three apples, a dozen bananas, and one cantaloupe for $2.36.
%Bob buys a dozen apples and two cantaloupes for $5.26. Carol buys two
%bananas and three cantaloupes for $2.77. How much do single pieces of
%each fruit cost?
A=[3 12 1;12 0 2;0 2 3]
b=[2.36 5.26 2.77]'

x=A\b

%% Forward Slash
x=b'/A'


%% Inconsistent Singular matrix
%Assume now that Carol buys six apples and
%one cantaloupe for $2.77.

A=[3 12 1;12 0 2;6 0 1]
b=[2.36 5.26 2.77]'

x=A\b

%% Consistent matrix
% The source of the difficulty is that the new information about Carol’s purchase
% is inconsistent with the earlier information about Alice’s and Bob’s purchases. We
% have said that Carol bought exactly half as much fruit as Bob. But she paid 2.77
% when half of Bob’s payment would have been only 2.63. The third row of A is equal
% to one-half of the second row, but b(3) is not equal to one-half of b(2). For this
% particular matrix A and vector b, the solution to the linear system of equations
% Ax = b does not exist.
% What if we make Carol’s purchase price consistent with Bob’s? We leave A
% unchanged and revise b with

b(3)=2.63

%% one partiluar solution
U=rref([A b])

%% Let x(3) = 0,solve Ax=b
Achange = [3 12 0;12 0 0;6 0 0]
Xp=[2.63/6 (2.36-2.63/2)/12 0]'
A*Xp

%% calc Null space of A
z=null(A,'r')
A*z

%% general solution of Ax=0
t=rand
Xn=t*z

%% Ax=b的通解等于Ax=b的特解 + Ax=0的通解
x=Xn+Xp
A*x

补充部分:

这里我也把求解Ax=0和Ax=b的全部过程写了下来:

解Ax=0:

在求解Ax=0时,我们要注意Ax=0的解等于高斯消元后的矩阵Rx=0的解,也就是说,高斯消元(行与行之间的消元)不改变A的零空间。

解Ax=b: 

 (全文完) 

作者 --- 松下J27

格言摘抄:

《长歌行》

青青园中葵,朝露待日晞。
阳春布德泽,万物生光辉。
常恐秋节至,焜黄华叶衰。
百川东到海,何时复西归?
少壮不努力,老大徒伤悲!

鸣谢(参考文献):

1,《Experiments with MATLAB》 - Cleve Moler

 (*配图与本文无关*) 

版权声明:所有的笔记,可能来自很多不同的网站和说明,在此没法一一列出,如有侵权,请告知,立即删除。欢迎大家转载,但是,如果有人引用或者COPY我的文章,必须在你的文章中注明你所使用的图片或者文字来自于我的文章,否则,侵权必究。 ----松下J27

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐