第五讲 矩阵与线性方程模型

第五讲矩阵与线性方程模型本章主要涉及线性代数中向量和矩阵的运算、线性方程组的解法等知识，介绍线性代数再经济方面的应用，以及MATLAB有关线性代数运算的命令。6.1引例：某城镇有三个重要产业，一个煤矿，一个发电厂和一条地方铁路。开采一元钱的煤，煤矿要支付0.25元的电费和0.35元的运输费；生产一元钱的电，发电厂要支付0.40元的煤费，0.05元的电费及0.10元的运输费；提供一元钱的运输服务，铁路要支付0.45元的煤费的0.10元的电费和0.10元的运输费，在某一周内煤矿接到外地金额50000元定货，发电厂接到外地金额25000元定货，外界对地方铁路需求为30000元.问:（1）三个企业间一周内总产值多少才能满足自身及外界需求？（2）三个企业间相互支付多少金额？三个企业各创造多少新价值？（3）如果煤矿需要增加总产值10000元，它对各个企业的产品或服务的完全需求分别将是多少？（4）假定三企业的外部需求仍是用于城镇的各种消费和积累，其中用于消费的产品价值分别为35000元、18000元和20000元，而假定三个企业的新创造价值又包括支付劳动报酬（工资等）和纯收入，其中支付劳动报酬分别为25488元、10146元和14258元，试分析各企业产品使用情况的比例关系；以及该星期系统的经济效益；（5）若在以后的三周内，企业外部需求的增长速度分别是15%、3%和12%；那么各企业的总产值将增长多少？投入产出综合平衡分析：在一个国家或区域的经济系统中，各部门（或企业）既有投入又有产出。生产的产品满足系统内部各部门和系统外的需求，同时也消耗系统内各部门内的产品。应如何组织生产呢？俄裔美国经济学家W.Leontief于20世纪30年代首先提出并成功地建立了研究国民经济投入产出地数学模型，他数次主持制定了美国的国民经济投入产出列表，且由此对国民经济各部门的结构和各种比例关系进行了定量分析。这一方法即投入产出法，投入产出法以其重要的应用价值迅速为世界各国经济学界和决策部门所采纳。W.Leontief因此于1973年获得了Nobel经济学奖。数学模型：假设有n个经济部门，xi为部门i的总产出；aij为部门j单位产品对部门i产品的消耗；yi为外部对部门i的需求；zj为部门j新创造的价值。则（6.1）称为分配平衡方程组，每一个等式以价值形式说明了对每一部门：中间产品（作为系统内个部门的消耗）+最终产品（外部需求）=总产品（6.2）称为消耗平衡方程组，说明对每一部门：对系统内各企业产品的消耗+新创价值=总价值11,1,,(6.1),1,,(6.2)niijjijnjjijjixaxyinxxazjn令（6.1）化为X=AX+Y(6.3)令C=E-A,E—单位矩阵，(6.3)化为CX=Z(6.4)经济学上称：A—直接消耗矩阵，C—Leontief矩阵，令B为投入产出矩阵，D为总投入向量，F为新创价值向量。111()(,,)'(,,)'(,,)'ijnnnAaXxxYyyZzz1,1,,1(6.5)nxBCDBFXDx6.2线性代数基本知识复习线性方程组就是n元一次方程。实际的问题，如大型的土建结构、机械结构、大型的输电网络、管道网络、投入产出分析都可以归结为线性方程组，有些数值计算的方法也会导致线性方程组的求解，如数据拟合、非线性方程组和偏微分方程组数值解问题等等。n往往达到几百、几千甚至上万。6.2.1线性方程组n个未知量m个方程的线性方程组的一般形式为1111111(6.6)nnmmnnmaxaxbaxaxb令则得矩阵形式(6.7)若b=0，即Ax=0(6.8)则称方程组为齐次的。方程组(6.7)可能有唯一解、有无穷多解，也可能无解，这取决于系数矩阵A及增广矩阵(A,b)的秩。若秩(A)=秩(A,b)=n，存在唯一解。其理论上可用Cramer法则求出，但计算量太大。若秩(A)=秩(A,b)n，存在无穷多解，其通解可表示为对应齐次方程组(6.8)的一个基础解系与(6.7)式的一个特解的叠加；若秩(A)≠秩(A,b)，则无解，这时一般寻求最小二乘近似解，即求x使向量Ax-b长度最小。111111,,nmmnnmxaabAxbaaxbAxb6.2.2逆矩阵方阵A称为可逆的，如果存在方阵B，使AB=BA=E这里E表示单位阵。B称为A的逆矩阵，记。A可逆的充分必要条件是。且(6.9)这里A*为A的伴随矩阵。(6.7)可表示为由于这个公式涉及大量的行列式计算，数值计算不采用。1BA0A11*AAA1xAb6.2.3特征值与特征向量对于方阵A，若存在数和非零向量x使称为A的一个特征值，x为A的一个对应于特征值的特征向量。特征值的计算归结为特征多项式的求根。对于n阶实数方阵，特征多项式在复数范围内总有n个根。对应于特征值的特征向量是齐次线性方程组的所有非零解。通常只需要求它的一组线性无关解。特征值和特征向量求解的数值方法是相当复杂的，适用性较广的是正交三角分解系列算法。(6.10)Axx()(6.11)0AEx6.3线性代数运算的MATLAB命令MATLAB是矩阵化程序设计语言，所以处理矩阵和向量运算特别方便。关于矩阵和向量的一些基本运算命令已在前面有所介绍，常用的命令和函数还有zeros生成0矩阵eig特征值、特征向量ones生成1矩阵diag对角矩阵eye生成单位矩阵trace方阵的迹linspace生成等距行向量rank矩阵的秩rand生成随机矩阵rref行最简形det方阵的行列式orth正交规范inv方阵的逆null求基础解系norm范数jordanJordan分解cond方阵的条件数6.3.1生成特殊矩阵zeros(m,n)m行n列0矩阵；ones(m,n)m行n列1矩阵；rand(m,n)m行n列[0，1]上均匀分布随机矩阵；eye(n)n阶单位矩阵；diag(A)A的对角线构成的向量（A为矩阵）；diag(X)X的元素构成的对角矩阵（X为向量）；linspace(x1,x2,n)x1与x2间的n维等距向量，即将[x1,x2]n-1等分。例：ones(4,4)ans=1111111111111111eye(4)ans=1000010000100001linspace(1,3,4)ans=1.00001.66672.33333.0000rand(2,4)%由于随机性，每次结果不同ans=0.95010.60680.89130.45650.23110.48600.76210.01856.3.2行列式和逆矩阵det(A)返回方阵A的行列式inv(A)返回A的逆矩阵例：A=[4,6;5,7];det(A),inv(A)ans=-2ans=-3.50003.00002.5000-2.00006.3.3矩阵除法左除法A\B解矩阵方程AX=B右除法B/A解矩阵方程XA=B（1）当A为方阵，其结果与inv(A)*B基本一致；（2）当A不为方阵，除法将分三种情况自动检测：若为超定方程组（既无解）除法将给出最小二乘意义上的近似解，即使向量AX-B的长度最小；若为不定方程组（即无穷多解），除法将给出一个具有最多零元素的特解（不是通解）；若为唯一解，除法将给出这个解。用户对结果应有一个正确的认识。例：解下列方程组11(4212324213324(2214242xyxyxyzxyzxyxyxyxyxy（）定解方程组）（）（不定方程组）（）超定方程组）（）（奇异方程组）解：A=[11;1-1];B=[1;4];x=A\Bx=2.5000-1.5000求得唯一解。A=[121;3-21];B=[1;4];x=A\Bx=1.2500-0.12500仅求得一个特解。A=[12;3-2;1-1];B=[1;4;2];x=A\Bx=1.2838-0.1757求得一最小二乘近似解。A=[12;24];B=[1;2];x=A\BWarning:Matrixissingulartoworkingprecision.(TypewarningoffMATLAB:singularMatrixtosuppressthiswarning.)x=InfInf可见，不能直接求解。A=[12;24;00];B=[1;2;0];x=A\B%增加0x+0y=0，使A不为方阵Warning:Rankdeficient,rank=1tol=2.9790e-015.x=00.5000仍可求一特解。例：求线性方程组的通解解：在有无穷多解的情况可用三种方法求得通解。12341234123411221xxxxxxxxxxxx方法一：用rref化为行最简形以后求解。clear;a=[1-11-1;-111-1;2-2-11];b=[1;1;-1];[rank(a),rank([a,b])]ans=22％秩相等且小于４，说明有无穷多解rref([a,b])ans=1-1000001-1100000即通解为：小x1=x2，x3=x4+1(x2，x4自由)方法二：先用除法求出一个特解，再用null求得齐次组的基础解系。clear;a=[1-11-1;-111-1;2-2-11];b=[1;1;-1];x0=a\b;x=null(a)Warning:Rankdeficient,rank=2tol=2.1756e-015.x=-0.70710-0.70710-0.00000.7071-0.00000.7071通解为k1*x(:,1)+k2*x(:,2)+x0方法三：使用solve求解。（见第８章）6.3.4特征值和特征向量[V,D]=eig(A)返回方阵A的特征值和特征向量。其中D为特征值构成的对角阵，每个特征值对应的V的为属于该特征值的一个特征向量，每个特征向量都是单位向量，并且属于同一特征值的线性无关特征向量已正交化。eig(A)返回方阵A的特征值构成的列向量。例：A=[123;234;245];[V,D]=eig(A),t=eig(A)V=-0.3957-0.2167+0.5832i-0.2167-0.5832i-0.57650.63130.6313-0.7149-0.3914-0.2471i-0.3914+0.2471iD=9.3329000-0.1665+0.2818i000-0.1665-0.2818it=9.3329-0.1665+0.2818i-0.1665-0.2818i6.4引例分析6.4.1数学模型记：x1－本周内煤矿总产值；x2－电厂总产值；x3－铁路总产值，则若记则12311232123300.400.45500000.250.050.10250000.350.150.1030000xxxxxxxxxxxx12350000,25000,3000000.400.450.250.050.1030.100.10xXxYxAAXYX此既说明了对每一企业：中间产品（作为系统内各企业的消耗）+最终产品（外部需求）=总产品称为分配平衡方程组。另一方面，若设z1，z2和z3（元）分别为煤矿、电厂和铁路在这个星期的创新价值，那么应有此说明对每一企业：对系统内各企业产品的消耗+创新价值+总产值称为消耗平衡方程组。11111222223333300.250.350.400.050.10(6.13)0.450.100.10xxxzxxxxzxxxxzx6.4.2问题的分析与解决Ⅰ.各需求量和新创价值在经济学上A称之为直接消耗矩阵，X称之为产出向量，Y称之为最后需求向量；A的元素称之为直接消耗系数。(6.12)又可以写为(E-A)X=Y(6.