MATLAB数学建模与仿真（第2版·微课视频版）程序代码 第9章-Matlab在线性代数中的应用

【例9-1】求下列矩阵的逆。>>a=[12;21];>>inv(a)ans=-0.33330.66670.6667-0.3333>>a^-1ans=-0.33330.66670.6667-0.3333【例9-2】求下列矩阵的逆。>>a=[21-3-1;3107;-124-2;10-15];>>inv(a)ans=-0.04710.5882-0.2706-0.94120.3882-0.35290.48240.7647-0.22350.2941-0.0353-0.4706-0.0353-0.05880.04710.2941【例9-3】计算行列式值，相应的MATLAB代码为。>>D=[31-12;-513-4;201-1;1-53-3];>>det(D)ans=40.0000【例9-4】计算行列式的值，相应的MATLAB代码为。>>symsabcd;>>D=[abcd;aa+ba+b+ca+b+c+d;a2*a+b3*a+2*b+c4*a+3*b+2*c+d;…>>a3*a+b6*a+3*b+c10*a+6*b+3*c+d];>>det(D)ans=a^4【例9-5】求行列式A=的值。>>A=[1,2;3,4]；>>A=[1,2;3,4];det(A)ans=-2【例9-6】求向量X=与向量Y=的点乘。>>X=[-102];>>Y=[-2-11];>>Z=dot(X,Y)Z=4还可用另一种算法：>>sum(X.*Y)ans=4【例9-7】利用函数rand和函数round构造一个5×5的随机正整数矩阵A和B。（1）计算A＋B，A－B和6A。>>A=round(rand(5)*10)A=121783104075109841058971081072>>B=round(rand(5)*10)B=78453074473386701087210281>>C=A+BC=81051211317841481317141110151616911812153>>D=A-BD=-6-6-325330-402712-310-5025988-11>>E=6*AE=61264248186024042306054482460304854426048604212（2）计算，和。>>(A*B)'ans=18406610110295106217229236929820421021213213520723720746114168146172>>(B)'*(A)'ans=18406610110295106217229236929820421021213213520723720746114168146172>>(A*B)^100ans=1.0e+287*0.28220.73950.68200.75950.55420.36220.94920.87550.97490.71130.62171.62921.50261.67321.22080.65241.70951.57671.75571.28100.65441.71481.58151.76111.2850（3）计算行列式，和。>>det(A)ans=4.0620e+03>>det(B)ans=2.1550e+04>>det(A*B)ans=8.7536e+07（4）若矩阵A和B可逆，计算和。>>inv(A)ans=0.33550.0487-0.4424-0.36040.63340.44290.0768-0.3183-0.59820.6901-0.9121-0.06650.78511.0369-1.31830.3904-0.0768-0.1817-0.40180.4766-0.25530.05020.19570.4165-0.5039>>inv(B)ans=0.1283-0.07470.0402-0.0624-0.01840.04460.0659-0.09380.0463-0.0310-0.0508-0.14150.14940.0949-0.0926-0.00150.0269-0.0468-0.00360.1507-0.01480.14250.0352-0.0985-0.0022【例9-8】向量a=(1,2,3),b=(2,3,4),c=(5,2,1),求a·(b×c)的混合积。>>a=[123]a=123>>b=[234]b=234>>c=[521]c=521>>v=dot(a,cross(b,c))v=-2【例9-9】求向量组a1=(1,-2,2,3),a2=(-2,4,-1,3),a3=(-1,2,0,3),a4=(0,6,2,3)的秩，并判断其线性相关性。>>A=[1-223;-24-13;-1203;0623];>>rank(A)ans=3由于秩为3，小于向量的个数，因此向量组线性相关。【例9-10】设M=，求矩阵M的秩。>>M=[32-1-3-2;2-131-3;705-1-8];>>rank(M) ans=2 【例9-11】已知矩阵M=QUOTE32-1-32-13170t-1秩为2，求常数t的值。解：题中的二阶子式不等于0.由于矩阵的秩为2，因此其三阶子式应该等于0。>>symst>>M=[32-1-3;2-131;70t-1];>>det(M(1:3,1:3))ans= 35-7*t 当t=5时，有一个三阶子式等于0，但是否所有的三阶子式都为0呢？输入：>>M=[32-1-3;2-131;705-1];>>rank(M)ans=2说明此时矩阵的秩为2。【例9-12】求下列矩阵列向量组a=的一个最大无关组。解：编写M文件ex1.m如下：formatrat；a=[1-,-102;-2426-6;2-1023;33334];b=rref(a)b=101/3016/3012/30-1/90001-1/300000由此可以得到5个列向量，并且易得向量组中的一个最大无关组。【例9-13】求向量组a=(1,-2,2,3),b=(-2,4,-1,3),c=(-1,2,0,3),d=(0,6,2,3),e=(2,-6,3,4)的一个最大无关组。>>a=[1-223]';>>b=[-24-13]';>>c=[-1203]';>>d=[0623]';>>e=[2-634]';>>A=[abcde]A=1-2-102-2426-62-102333334>>formatrat%以有理格式输出>>B=rref(A)%求A的行最简形B=101/3016/9012/30-1/90001-1/300000从B中可以得到：向量组a、b、d是其中一个最大无关组。【例9-14】求向量组a=(1,-1,2,4)，b=(0,3,1,2)，c=(3,0,7,14)，d=(1,-1,2,0)，E=(2,1,5,0)的最大无关组，并将其他向量用最大无关组线性表示。>>a=[1-124;0312;30714;1-120;2150];>>b=transpose(a);>>rref(b)ans=1.000003.00000-0.500001.00001.000001.00000001.00002.500000000在行最简形中有3个非零行，因此向量的秩等于3.非零行的首元素分别位于第一，二，四列，因此a,b,d是向量组的一个最大无关组。第三列的前俩个元素分别是3,1。于是c=3a+b。第5列的前3个元素分别是-0.5,1,2.5。于是e=-0.5a+b+2.5d。【例9-15】已知向量组，，，，，求出它的最大无关组，并用该最大无关组来线性表示其它向量。>>formatrat>>A=[31290;41338;0002-2;826121;321210]A=31290413380002-2826121321210>>B=rref(A)B=1010201-1030001-10000000000,,是向量组的一个最大无关组。且有=-，=2+3–。【例9-16】求下列向量组的秩及一个最大线性无关组，并将其余向量用最大线性无关组表示，。分析：容易发现用定义的形式很难求秩和最大线性无关组，为此我们从方程组和矩阵之间的关系以及方程组和向量组之间的关系可以得到，向量组的秩及其最大线性无关组应该与其对应的矩阵的秩以及矩阵的最高阶非零子式之间有某种关系，为此我们给出：定理：矩阵的秩等于其行向量组的秩，也等于其列向量组的秩。略证：设A的秩为r，则在A中存在r阶子式，从而所在的r线性无关，又A的所有的r+1阶子式，因此中的任意个列向量都线性相关，因此所在的列是A列向量组的最大线性无关组，所以列向量组的秩等于r。类似可证矩阵A行向量组的秩等于r。同时从证明的过程可以发现：若是矩阵A的一个最高阶非零子式，则所在的r列即是A的列向量组的一个最大线性无关组；同时所在的r行即是A的行向量组的一个最大线性无关组。我们现在求解上面的问题，把上面的4个向量看成某矩阵A的4列进行求解。解：所以，是的一个最大线性无关组。（当然易见亦是的一个最大线性无关组）。为了把用线性表示，把A再变成行最简形矩阵易见。初等变换前后列向量组之间的线性表示形式是保持不变的，同时可以验证上面的线性表示的结果是正确的。【例9-17】求方程组的一个特解。>>A=[11-3-1;3-1-34;15-9-8];>>b=[140]';>>[Q,R]=qr(A)Q=-0.30150.1421-0.9428-0.9045-0.35530.2357-0.30150.92390.2357R=-3.3166-0.90456.3317-0.904505.1168-7.6752-8.954400-0.0000-0.0000>>X=R\(Q\b)Warning:Rankdeficient,rank=2tol=8.8373e-015.X=00-0.53330.6000说明：这3种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。【例9-19】求解方程组的通解。>>A=[1221;21-2-2;1-1-4-3];>>formatrat>>B=null(A,'r')B=25/3-2-4/31001写出通解：>>symsk1k2>>x=k1*B(:,1)+k2*B(:,2)%写出方程组的通解>>symsk1k2>>x=k1*B(:,1)+k2*B(:,2)x=2*k1+(5*k2)/3-2*k1-(4*k2)/3k1k2>>pretty(x)%让通解表达式更加精美/5k2\|2k1+----||3||||4k2||-2k1-----||3||||k1|||\k2/【例9-20】求解方程组的解。在MATLAB编辑器中建立M文件：LX0601.mA=[1-23-1;3-15-3;212-2];b=[123]';B=[Ab];n=4;R_A=rank(A);R_B=rank(B);formatratifR_A==R_B&R_A==n%判断有唯一解X=A\belseifR_A==R_B&R_A<n%判断有无穷解X=A\b%求特解C=null(A,'r')%求AX=0的基础解系elseX='equationnosolve'%判断无解End运行后结果显示：X=equitionnosolve说明该方程组无解。【例9-20】求解方程组的通解。方法一：在MATLAB编辑器中建立M文件：LX0602.mA=[11-3-1;3-1-34;15-9-8];b=[140]';B=[Ab];n=4;R_A=rank(A);R_B=rank(B);formatratifR_A==R_B&R_A==nX=A\belseifR_A==R_B&R_A<nX=A\bC=null(A,'r')elseX='equationhasnosolves'end运行结果显示为：>>LX0602.m警告:秩亏，秩=2，tol=8.837264e-15。>InLX0602.m(line11)X=00-8/153/5C=3/2-3/43/27/41001所以原方程组的通解为>>symsk1k2>>X=k1*C(:,1)+k2*C(:,2)+XX=(3*k1)/2-(3*k2)/4(3*k1)/2+(7*k2)/4k1-8/15k2+3/5>>pretty(X)/3k13k2\|---------||24||||3k17k2||----+----||24||||8||k1---||15||||3||k2+-|\5/方法二：在MATLAB编辑器中建立M文件：LX0603.mA=[11-3-1;3-1-34;15-9-8];b=[140]';B=[Ab];C=rref(B)%求增广矩阵的行最简形，可得最简同解方程组。运行结果：C=10-3/23/45/401-3/2-7/4-1/400000对应齐次方程组的基础解系为：非齐次方程组的特解为：所以原方程组的通解为：【例9-21】求行列式矩阵A=的特征向量矩阵，特征量矩阵。>>A=[123;456;789];>>[x,y]=eig(A)x=-0.2320-0.78580.4082-0.5253-0.0868-0.8165-0.81870.61230.4082y=16.1168000-1.1168000-0.0000其中x为特征向量矩阵，y为特征值矩阵。【例9-22】求矩阵的特征值和特征向量。相应的MATLAB代码和计算结果如下：>>A=[3-1;-13]A=3-1-13>>eig(A)%A的特征值ans=24>>[X,D]=eig(A)%D的对角线元素是特征值，X是矩阵X=-985/1393-985/1393-985/1393985/1393D=2004【例9-23】求下列矩阵的特征值和特征向量，并判断其正定性。（1）>>A=[123;256;3625]A=1232563625>>formatrat>>[V,D]=eig(A)V=160/171445/13571377/10567-751/21351596/1781417/1541-301/10736-712/2381909/953D=25/1580003767/10100003145/116即特征值25/158对应特征向量（160/171,-751/2135,-301/10736）T；特征值3767/1010对应特征向量（445/1375,1596/1781,-712/2381）T；特征值3145/116对应特征向量（1377/10567,417/1541,909/953）T。因为A的特征值均为正数，所以A正定。（2）>>B=[-2031;3-10-6;1-6-22]B=-20313-10-61-6-22>>[V,D]=eig(B)V=-357/9374822/5323500/27031060/2647-19/10193681/40187996/9595699/1652-1609/4524D=-20323/802000-7348/375000-544/77特征值与特征向量的对应关系如上，因为B的特征向量均为负数，所以B负定。【例9-24】求矩阵A=[61219；-9-20-33；4915]的特征值,具体代码序列如下。>>A=[61219;-9-20-33;4915];>>d=eig(A)d=-1.00001.00001.0000【例9-25】将矩阵A=正交化。>>A=[400;031;013];>>P=orth(A)P=010-985/13930-985/1393-985/13930985/1393>>formatshort>>PP=01.00000-0.70710-0.7071-0.707100.7071>>Q=P'*PQ=1.000000.000001.000000.000001.0000【例9-26】将a1=，a2=，a3=向量组正交化。>>a1=[1;0;-1;0]>>a1=10-10>>a2=[1;-1;0;1];>>a3=[-1;1;1;0];>>Q=orth([a1a2a3])Q=-0.69280.0587-0.42800.50460.4078-0.76090.4589-0.6730-0.0563-0.2339-0.6143-0.4843>>Q'*Qans=1.00000-0.000001.00000.0000-0.00000.00001.0000Q就是正交化后的矩阵，orth()是正交化函数【例9-27】求矩阵A=的schur分解。>>A=[400;031;013];>>[U,T]=schur(A)U=001.0000-0.70710.707100.70710.70710T=200040004这里U就是所求的正交矩阵P，T就是对角矩阵Λ。>>[V,D]=eig(A)V=001.0000-0.70710.707100.70710.70710D=200040004这里V就是所求的正交矩阵P，D就是对角矩阵Λ。说明：对于实对称矩阵，用eig和schur分解效果一样。【例9-30】用一个正交变换X=PY，把二次型化成标准形。解：先写出二次型的实对称矩阵在MATLAB编辑器中建立M文件：LX0603.mA=[011-1;10-11;1-101;-1110];[P,D]=schur(A)运行结果：>>LX0603P=-1/2390/1351780/989780/36911/2-390/1351780/3691780/9891/2-390/1351780/1351-780/1351-1/2-1170/135100D=-3000010000100001继续输入：>>symsy1y2y3y4>>y=[y1;y2;y3;y4];%这里不用行向量的转置表示，以免出现复数。>>X=vpa(P,2)*y%vpa表示可变精度计算，这里取2位精度。X=0.28867513459408655762672424316406*y2-0.5*y1+0.78867513465229421854019165039062*y3+0.21132486540591344237327575683594*y40.5*y1-0.28867513459408655762672424316406*y2+0.21132486540591344237327575683594*y3+0.78867513465229421854019165039062*y40.5*y1-0.28867513459408655762672424316406*y2+0.57735026918817311525344848632812*y3-0.57735026918817311525344848632812*y4-0.5*y1-0.86602540384046733379364013671875*y2>>f=[y1y2y3y4]*D*yf=-3*y1^2+y2^2+y3^2+y4^2【例9-29】利用U=[U,L]=schur(A)计算A=[1,4,2;5,6,9;4,1,8]的schur分解，具体代码如下：>>A=[1,4,2;5,6,9;4,1,8];>>[U,L]=schur(A)U=0.34940.89290.28380.8242-0.1489-0.54640.4456-0.42480.7880L=12.9859-0.38997.27070-0.4658-3.9981002.4799【例9-30】求一个正交变换X=PY，把二次型化成标准形。先写出二次型的实对称矩阵：在MATLAB编辑器中建立M文件如下：A=[011-1;10-11;1-101;-1110];[P,D]=schur(A)symsy1y2y3y4y=[y1;y2;y3;y4];X=vpa(P,2)*y%vpa表示可变精度计算，这里取2位精度f=[y1y2y3y4]*D*y运行后结果显示如下：P=-1/2390/1351780/989780/36911/2-390/1351780/3691780/9891/2-390/1351780/1351-780/1351-1/2-1170/135100D=-3000010000100001X=0.28867513459408655762672424316406*y2-0.5*y1+0.78867513465229421854019165039062*y3+0.21132486540591344237327575683594*y40.5*y1-0.28867513459408655762672424316406*y2+0.21132486540591344237327575683594*y3+0.78867513465229421854019165039062*y40.5*y1-0.28867513459408655762672424316406*y2+0.57735026918817311525344848632812*y3-0.57735026918817311525344848632812*y4-0.5*y1-0.86602540384046733379364013671875*y2f=-3*y1^2+y2^2+y3^2+y4^2【例9-31】求一个正交变换，把二次型化为标准型,二次型的矩阵为：A=。代码如下：>>A=[0,1,1,-1;1,0,-1,1;1,-1,0,1;-1,1,1,0];>>[P,D]=eig(A)P=-1/2390/1351780/989780/36911/2-390/1351780/3691780/9891/2-390/1351780/1351-780/1351-1/2-1170/135100D=-3000010000100001P就是所求的正交矩阵，使得PTAP=D，令X=PY，其中X=[x1,…,x4]T，Y=[y1,…,y4]T化简后的二次型为。上面求得的正交矩阵P是数值解，下面我们求正交矩阵的精确解。>>a=sym([011-1;10-11;1-101;-1110]);>>[v,d]=eig(a)求得：v=[1,1,1,-1][-1,1,0,0][-1,0,1,0][1,0,0,1]d=[-3,0,0,0][0,1,0,0][0,0,1,0][0,0,0,1]即求得矩阵的特征值为1、1、1、-3，对应的特征向量分别是矩阵v的第1、2、3、4列。再把对应于特征值1的3个特征向量正交化、单位化，我们就容易求出正交矩阵。【例9-32】求线性方程组的一个特解和一个基础解系。>>formatrat>>a=[11100;11-1-1-2;220-1-2;55-3-4-8];>>b=[0;1;1;4];>>x=rank(a);>>y=rank([a,b]);>>n=size(a,2);>>ifx==y&x==nfprintf('方程组有唯一解\n')x=a\belseifx==y&x<n fprintf('方程组有无穷多个解\n')x=a\b%求非齐次方程组的特解xt=null(a,'r')%求齐次方程的基础解系end结果如下：方程组有无穷多个解警告:秩亏，秩=2，tol=9.420555e-15。x=0-1/166331309922504500-1/2xt=-11/211000-1/2-1010001【例9-33】用MATLAB解下列线性方程组。>>A=[24-6;153;132];>>b=[-4;10;5];>>x=inv(A)*bx=-321>>B=[341-62;4503;113825];>>c=[-41;100;50];>>x=inv(B)*cx=-1835/208611/2361964/1009【例9-34】求解下面的方程组：分析：对于线性方程组求解，常用线性代数的方法，把方程组转化为矩阵进行计算。MATLAB的表达形式及结果如下：>>a=[816;357;492];%建立系数矩阵>>b=[7.5;4;12];%建立常数项矩阵>>x=a\b%求方程组的解计算结果：x=931/720323/360-449/720【例9-35】本金P以每年n次，每次i%的增值率（n与i的乘积为每年增值额的百分比）增加，当增加到r×P时所花费的时间T为：（利用复利计息公式可得到下式）。（）MATLAB的表达形式及结果如下：>>r=2;i=0.5;n=12;%变量赋值>>T=log(r)/(n*log(1+0.01*i))T=3347/289计算结果显示为：T=3347/289即所花费的时间为T=3347/289年。分析：上面的问题是一个利用公式直接进行赋值计算问题，实际中若变量在某个范围变化取很多值时，使用MATLAB，将倍感方便，轻松得到结果，其绘图功能还能将结果轻松的显示出来，变量之间的变化规律将一目了然。若r在[1,9]变化，i在[0.5,3.5]变化，将MATLAB的表达式作如下改动，结果如图9-1。>>r=1:0.5:9;>>i=0.5:0.5:3.5;>>n=12;>>p=1./(n*log(1+0.01*i));>>T=log(r')*p;>>plot(r,T)>>xlabel('r')%给x轴加标题>>ylabel('T')%给y轴加标题>>q=ones(1,length(i));>>text(7*q-0.2,[T(14,1:5)+0.5,T(14,6)-0.1,T(14,7)-0.9],num2str(i'))图9-1T-r曲线从图9-1中既可以看到T随r的变化规律，而且还能看到i的不同取值对T-r曲线的影响（图中的7条曲线分别代表i的不同取值）。【例9-36】人口迁徙模型问题。设在一个大城市中的总人口是固定的。人口的分布则因居民在市区和郊区之间迁徙而变化。每年有6%的市区居民搬到郊区去住，而有2%的郊区居民搬到市区。假如开始时有30%的居民住在市区，70%的居民住在郊区，问十年后市区和郊区的居民人口比例是多少？30年、50年后又如何？这个问题可以用矩阵乘法来描述。把人口变量用市区和郊区两个分量表示，即其中xc为市区人口所占比例，xs为郊区人口所占比例，k表示年份的次序。在k=0的初始状态：。一年以后，市区人口为；郊区人口，用矩阵乘法来描述，可写成：此关系可以从初始时间到k年,扩展为,用下列MATLAB程序进行计算：>>A=[0.940.02;0.060.98];>>x0=[0.3;0.7];>>x1=A*x0;>>x10=A^10*x0;>>x30=A^30*x0;>>x50=A^50*x0;>>[x1,x10,x30,x50]ans=37/1251348/496131/1221135/452688/125394/54191/122242/323程序运行的结果转换为小数为：无限增加时间k，市区和郊区人口之比将趋向一组常数0.25/0.75。为了弄清为什么这个过程趋向于一个稳态值，我们改变一下坐标系统。在这个坐标系统中可以更清楚地看到乘以矩阵A的效果。选u1为稳态向量[0.25,0.75]T的任意一个倍数，令u1=[1,3]T和u2=[-1,1]T。可以看到，用A乘以这两个向量的结果不过是改变向量的长度，不影响其相角（方向）：初始向量x0可以写成这两个基向量u1和u2的线性组合；因此式中的第二项会随着k的增大趋向于零。如果只取小数点后两位，则只要k>27，这第二项就可以忽略不计而得到适当选择基向量可以使矩阵乘法结果等价于一个简单的实数乘子，避免相角项出现，使得问题简单化。这也是方阵求特征值的基本思想。这个应用问题实际上是所谓马尔可夫过程的一个类型。所得到的向量序列x1,x2,...,xk称为马尔可夫链。马尔可夫过程的特点是k时刻的系统状态xk完全可由其前一个时刻的状态xk-1所决定，与k-1时刻之前的系统状态无关。【例9-37】线性规划问题。线性规划问题即目标函数和约束条件均为线性函数的问题。其标准形式为：minC’xsub.ToAx=bx≥0其中C,b，，均为数值矩阵，。若目标函数为：maxC’x，则转换成：min–C’x。标准形式的线性规划问题简称为LP（LinearProgramming）问题。其它形式的线性规划问题经过适当的变换均可以化为此种标准形。线性规划问题虽然简单，但在工农业及其他生产部门中应用十分广泛。在MATLAB中，线性规划问题由linprog函数求解。函数：linprog%求解如下形式的线性规划问题：suchthat其中f,x,b,beq,lb,ub为向量，A,Aeq为矩阵。格式：x=linprog(f,A,b)x=linprog(f,A,b,Aeq,beq)x=linprog(f,A,b,Aeq,beq,lb,ub)x=linprog(f,A,b,Aeq,beq,lb,ub,x0)x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)[x,fval]=linprog(...)[x,fval,exitflag]=linprog(...)[x,fval,exitflag,output]=linprog(...)[x,fval,exitflag,output,lambda]=linprog(...)说明：x=linprog(f,A,b)求解问题minf’*x，约束条件为A*x<=b。x=linprog(f,A,b,Aeq,beq)求解上面的问题，但增加等式约束，即Aeq*x=beq。若没有不等式存在，则令A=[]、b=[]。x=linprog(f,A,b,Aeq,beq,lb,ub)定义设计变量x的下界lb和上界ub，使得x始终在该范围内。若没有等式约束，令Aeq=[]、beq=[]。x=linprog(f,A,b,Aeq,beq,lb,ub,x0)设置初值为x0。该选项只适用于中型问题，默认时大型算法将忽略初值。x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)用options指定的优化参数进行最小化。[x,fval]=linprog(...)返回解x处的目标函数值fval。[x,fval,exitflag]=linprog(...)返回exitflag值，描述函数计算的退出条件。[x,fval,exitflag,output]=linprog(...)返回包含优化信息的输出变量output。[x,fval,exitflag,output,lambda]=linprog(...)将x处的Lagrange乘子返回到lambda参数。exitflag参数描述退出条件：·>0表示目标函数收敛于解x处；·=0表示已经达到函数评价或迭代的最大次数；·<0表示目标函数不收敛。output参数包含下列优化信息：·output.iterations迭代次数；·output.cgiterationsPCG迭代次数(只适用于大型规划问题)；·output.algorithm所采用的算法。lambda参数是解x处的Lagrange乘子。它有以下一些属性：·lambda.lower—lambda的下界；·lambda.upper—lambda的上界；·lambda.ineqlin—lambda的线性不等式；·lambda.eqlin—lambda的线性等式。某厂生产甲乙两种产品，已知制成一吨产品甲需资源A3吨，资源B4m3；制成一吨产品乙需资源A2吨，资源B6m3；资源C7个单位。若一吨产品甲和乙的经济价值分别为7万元和5万元，三种资源的限制量分别为90吨、200m3和210个单位，试决定应生产这两种产品各多少吨才能使创造的总经济价值最高？解：令产品甲的数量为x1，产品乙的数量为x2。由题意可以建立下面的数学模型：sub.to该模型中要求目标函数最大化，需要按照MATLAB的要求进行转换，即目标函数为在MATLAB中实现：>>f=[-7;-5];>>A=[32;46;07];>>b=[90;200;210];>>lb=[0;0];>>[x,fval,exitflag,output,lambda]=linprog(f,A,b,[],[],lb)Optimalsolutionfound.x=1424fval=-218exitflag=1output=包含以下字段的struct:iterations:3constrviolation:*message:'Optimalsolutionfound.'algorithm:'dual-simplex'firstorderopt:*lambda=包含以下字段的struct:lower:[2×1double]upper:[2×1double]eqlin:[]ineqlin:[3×1double]由上可知，生产甲种产品14吨、乙种产品24吨可使创造的总经济价值最高为218万元。exitflag=1表示过程正常收敛于解x处。【例9-38】交通流的分析问题。某城市有两组单行道，构成了一个包含四个节点A,B,C,D的十字路口如图9-2所示。在交通繁忙时段的汽车从外部进出此十字路口的流量（每小时的车流数）标于图上。现要求计算每两个节点之间路段上的交通流量x1,x2,x3,x4。解：在每个节点上，进入和离开的车数应该相等，这就决定了四个流通的方程：节点A:x1+450＝x2+610节点B:x2+520＝x3+480节点C:x3+390＝x4+600节点D:x4+640＝x2+310将这组方程进行整理，写成矩阵形式：图9-2 图9-2 单行线交通流图其系数增广矩阵为：用消元法求其行阶梯形式，或者直接调用U0=rref([A,b]),可以得出其精简行阶梯形式为注意这个系数矩阵所代表的意义，它的左边四列从左至右依次为变量x1,x2,x3,x4的系数，第五列则是在等式右边的常数项。把第四列移到等式右边，可以按行列写恢复为方程，其结果为：x1=x4+330,x2=x4+170,x3=x4+2100＝0由于最后一行变为全零，这个精简行阶梯形式只有三行有效，也就是说四个方程中有一个是相依的，实际上只有三个有效方程。方程数比未知数的数目少，即没有给出足够的信息来唯一地确定x1,x2,x3,和x4。其原因也不难从物理上想象，题目给出的只是进入和离开这个十字路区的流量，如果有些车沿着这四方的单行道绕圈，那是不会影响总的输入输出流量的，但可以全面增加四条路上的流量。所以x4被称为自由变量，实际上它的取值也不能完全自由，因为规定了这些路段都是单行道，x1,x2,x3,和x4。都不能取负值。所以要准确了解这里的交通流情况，还应该在x1,x2,x3,和x4中，再检测一个变量。【例9-39】价格平衡模型问题。在Leontiff成为诺贝尔奖金获得者的历史中，线性代数曾起过重要的作用，我们来看看他的基本思路。假定一个国家或区域的经济可以分解为n个部门，这些部门都有生产产品或服务的独立功能。设单列n元向量x是这些n个部门的产出，组成在Rn空间的产出向量。先假定该社会是自给自足的经济，这是一个最简单的情况。因此各经济部门生产出的产品，完全被自己部门和其它部门所消费。Leontiff提出的第一个问题是，各生产部门的实际产出的价格p应该是多少，才能使各部门的收入和消耗相等，以维持持续的生产。Leontiff的输入输出模型中的一个基本假