MATLAB数学建模与仿真（第2版·微课视频版）程序代码汇 综合实例（1-10）

【例2-101】将多项式与三角函数结合。>>x=linspace(0,4*pi,10);>>y=sin(x);%在区间

[0,4*pi]

中沿正弦曲线生成10个等间距的点。>>p=polyfit(x,y,7);%使用

polyfit

将一个7次多项式与这些点拟合。>>x1=linspace(0,4*pi);>>y1=polyval(p,x1);>>figure运行以上程序代码后，得到如图2-4的图形。图2-4将多项式与三角函数结合的图形【例2-102】简单线性回归。>>x=1:50;y=-0.3*x+2*randn(1,50);p=polyfit(x,y,1);%将一个简单线性回归模型与一组离散二维数据点拟合。创建几个由%样本数据点(x,y)

组成的向量。对数据进行一次多项式拟合。f=polyval(p,x);plot(x,y,'o',x,f,'-')legend('data','linearfit')%计算在x中的点处拟合的多项式p。用这些数据绘制得到的%线性回归模型。运行以上代码后，得到如图2-5的图形。图2-5简单线性回归的图形【例2-103】使用中心化和缩放改善数值属性>>year=(1750:25:2000)';pop=1e6*[791856978105012621544165025326122817011560]';T=table(year,pop);%创建一个由1750-2000年的人口数据组成的表，并绘制数据点。plot(year,pop,'o')%生成图2-6>>[p,~,mu]=polyfit(T.year,T.pop,5);%使用带三个输入的polyfit拟合一个使用中心%化和缩放的5次多项式，这将改善问题的数值属性。>>f=polyval(p,year,[],mu);>>holdon>>plot(year,f)>>holdoff运行以上代码，得到如图2-6、2-7的图形。图2-6人口数据点图2-7中心化和缩放改善数值属性【例3-93】已知长方体的表面积为12QUOTEa2a2(a>0)，求体积最大的长方体的长、宽、高并计算出其体积。(1)根据题中要求可知，设长方体的长为x,宽为y,高为z；所求的目标函数为f（x,y,z）=xyz；限制条件为g（x,y,z）=2（xy+xz+yz）=12QUOTEa2a2，即QUOTEμμ(x，y，z)=2（xy+xz+yz）-12QUOTEa2a2=0；引入拉格朗日乘子QUOTEλλ，构造拉格朗日函数L(x，y，z)=f（x，y，z）+QUOTEλλ[2（xy+xz+yz）-12QUOTEa2a2]。(2)根据以上分析，输入代码。其中QUOTEλλ使用s表示。>>symsxyzsa;>>L=x*y*z+s*(2*y*z+2*z*x+2*x*y-12*a^2);%建立拉格朗日方程>>Lx=diff(L,'x')%对x求导Lx=y*z+s*(2*y+2*z)>>Ly=diff(L,'y')%对y求导Ly=x*z+s*(2*x+2*z)>>Lz=diff(L,'z')%对z求导Lz=x*y+s*(2*x+2*y)Ls=diff(L,'s')Ls=-12*a^2+2*x*y+2*x*z+2*y*z>>[sxyz]=solve(Lx,Ly,Lz,Ls)s=-(2^(1/2)*a)/4(2^(1/2)*a)/4x=2^(1/2)*a-2^(1/2)*ay=2^(1/2)*a-2^(1/2)*az=2^(1/2)*a-2^(1/2)*a（3）根据以上结果且x、y、z应都大于零，因此在表面积固定的情况下，长方体是存在最大体积的，即x=y=z=QUOTE2a2a时，长方体的体积最大。输入以下代码求出长方体体积。>>V=x.*y.*zV=2*2^(1/2)*a^3【例4-77】绘制饼图和复数向量的二维统计分析图代码：subplot(1,2,1);pie([2347,1827,2043,3025]);title('bingtu');legend('q1','q2','q3','q4');subplot(1,2,2);compass([7+2.9i,2-3i,-1.5-6i]);title('xiangliangtu');执行结果如图4-85所示。图4-85饼图和复数向量的二维统计分析图【例4-78】使用scatter3函数绘制三维立体图代码：x=rand(1,20);y=rand(1,20);z=x+y;figure;subplot(121);scatter3(x,y,z) %绘制三维散点图title('空心点');subplot(122);scatter3(x,y,z,'r','filled'); %绘制三维散点图title('实心点');执行结果如图4-86所示。图4-86三维立体图【例5-43】此综合示例将演示如何在Python代码中调用MATLAB的M文件脚本。配置环境Python版本要求：本文使用的MATLAB支持Python版本2.7、3.6和3.7。安装：将包含Python解释器的文件夹添加到您的路径（如果尚未在该路径中）。找到MATLAB文件夹的路径。启动MATLAB，并在命令行窗口中键入matlabroot。复制matlabroot所返回的路径。要安装引擎API，请选择以下选项之一（您可能需要管理员特权才能执行这些命令）：在Windows操作系统提示符下- cd"matlabroot\extern\engines\python"pythonsetup.pyinstall在macOS或Linux操作系统提示符下-cd"matlabroot/extern/engines/python"pythonsetup.pyinstall实现Python调用Matlab的M脚本文件。 MATLAB脚本文件必须与Python文件放在同一目录下。新建MATLAB脚本文件count.m，代码如下：1+2新建python文件test.py，代码如下：importtimeimportmatlab.engineeng=matlab.engine.start_matlab()#matlab.engine.start_matlab()启动一个新的MATLAB进程，并返回Python的一个变量，它是一个MatlabEngine对象，用于与MATLAB过程进行通信。eng.count(nargout=0)time.sleep(10)打开CMD，进入M脚本文件和Python文件所在的目录，输入以下命令：python.\test.py得到结果如图5-23所示：图5-23在Python调用MATLAB脚本文件【例6-9】创建一个Simulink模型HelloWorld。1.打开Matlab，新建一个Simulink模型，同时将模型另存为HelloWorld.slx。2.打开模块库浏览器。在Sources库中添加2个Constant模块，将其中一个旋转至输出端向上。在MathOperations库中添加1个Sum模块。在Sinks库中添加1个Display模块。3.连接信号线。如图6-69所示。图6-69系统框图4.按下快捷键Ctrl+B对模型进行编译。此时，由于还没有对代码生成进行相关配置，系统会弹出错误提示，要求对代码生成进行相关配置。由于程序仿真为非连续仿真，所以在生成C代码时，Simulink会给出三个提示及对策：需要将仿真解算器(Solver)步长改为固定步长(fixed-step)。把输出目标文件类型改为可变步长目标文件。通过向导完成代码生成的相关设置。如图6-70所示。图6-70错误提示5.点击第一项后的Open按钮，在弹出模型参数设置对话框中，将Solverselection栏中的Type改为fixed-step。如图6-71所示。图6-71设置Solverselection6.为了简化其他代码生成的相关设置，点击错误提示中的第三项后的Open按钮，利用SimulinkCoderQuickStart向导进行配置。如图6-72所示。图6-72SimulinkCoderQuickStart向导7.一直点击Next按钮，采用默认配置即可。配置完成后，点击Finish按钮关闭SimulinkCoderQuickStart向导。如图6-73所示。图6-73SimulinkCoderQuickStart向导配置完成8.回到Simulink模型界面，按下快捷键Ctrl+B进行C代码自动生成。由于之前的SimulinkCoderQuickStart向导会默认打开“代码生成报告”功能，因此在代码生成结束后，会自动弹出网页形式的报告。如图6-74所示。图6-74网页形式的报告9.点击报告左边栏中的HelloWorld.h可以打开生成代码的头文件。但是，在对应的HelloWorld.c中找不到加法运算的相关代码。原因在于在模型中只是将变量显示在了Simulink的Display模块中。对于程序来说，并没有真正的输出。没有输出的代码会被Simulink优化掉。10.将模型中的Constant和Display分别用in模块和Out模块替换。如图6-75所示。图6-75模块替换后的系统框图11.画框选中所有模块，按下快捷键Ctrl+G或通过MULTIPLE选项卡→CREATE栏→CreateSubsystem，将选中模块集生成一个子系统，以便最终生成一个叫GetSum的C代码函数。如图6-76所示。图6-76生成子系统后的系统框图12.在模型中添加一个DataStoreMemory模块，添加这个模块相当于在C语言中声明一个新的全局变量A。添加DataStoreRead模块替换In1模块，读取A的数据，作为GetSum子系统的输入。添加DataStoreWrite模块替换Out1模块，将GetSum子系统的输出，写入到A中。如图6-77所示。图6-77模块替换后的系统框图13.为了生成一个单独的函数，需要将GetSum子系统，定义为atomicunit（原子子系统）。在GetSum子系统上右键→BlockParameters(Subsystem)→选项卡Main→勾选Treatasatomicunit→选项卡CodeGeneration→Functionpackaging选择Reusablefunction。如图6-78所示。图6-78GetSum子系统参数设置14.按下快捷键Ctrl+B生成代码。在代码生成报告的左侧链接处，或在matlab文件夹中可以看到生成的代码文件HellowWorld.c。15.打开HelloWorld.c文件，可以看到HelloWorld_GetSum就是模型中GetSum子系统对应的C代码函数。如图6-79所示。图6-79HelloWorld.c中的函数HelloWorld_GetSum【例7-39】函数逼近:设计一个BP神经网络,g(x)=1+sin(kπ/2x)，其中,分别令

k=2,3,6进行仿真,通过调节参数得出信号频率和隐含层节点之间,隐含层节点与函数逼近能力之间的关系。解：设频率参数

k=2,绘制要毕竟的非线性函数的目标曲线:k=2;p=[-1:.05:8];t=1+sin(k*pi/2*p);plot(p,t,'-');title('Non-linearFunctiontobeApproached');xlabel=('Time');ylabl=('Non-linearFunction');运行后得到目标曲线如下图所示:使用newff函数建立BP神经网络n=5;%设定隐含层神经元数目为5net=newff(minmax(p),[n,1],{'tansig','purelin'},'trainlm');%选择隐含层和输出层神经元传递函数为tansig和purelin,%网络训练算法采用Levenberg-Marquardt算法trainlm.%对于初始神经网络,可使用sim()函数观察网络输出y1=sim(net,p);figure;plot(p,t,'-',p,y1,':')title('RawOutput');xlabel('Time');ylabel('SimulationOutput----OriginalFunction')程序运行后所得到的神经网络输出曲线和原函数的对比图:使用newff函数建立神经网络时权值和阈值的初始化是随机的,故网络输出的结构很差,达不到逼近函数的目的。下面设置网络训练参数,应用train函数对网络进行训练:net.trainParam.Epochs==200;%设定网络训练步长为200net.trainParam.goal=0.2;%设定网络训练精度为0.2net=train(net,p,t);%开始执行网络训练下面对训练出的神经网络进行仿真：y2=sim(net,p);figure;plot(p,t,'-',p,y1,':',p,y2,'--')title('OutputoftheTrainedNeuralNetwork');clearxlabelxlabel('Time')ylabel('SimulationOutput')下图为训练所得的神经网络的输出结果:不难看出,经过训练后的BP神经网络对非线性函数的逼近效果有显著提升。实际上,改变非线性函数的频率和BP神经网络隐含层神经元的数目,也会对神经网络逼近函数的效果有一定影响.一般而言,隐含层神经元数目越多,BP网络逼近非线性函数的能力越强。【例8-38】实验一：生日蛋糕问题。一个数学家即将要迎来他九十岁生日，有很多的学生要来为他祝寿，所以要定做一个特大蛋糕。为了纪念他提出的一项重要成果——口腔医学的悬链线模型，他的弟子要求蛋糕店的老板将蛋糕边缘圆盘半径做成下列悬链线函数：r=2-(exp(2h)+exp(-2h))/5，0<h<1(单位m)由于蛋糕店从来没有做过这么大的蛋糕，蛋糕店的老板必须要计算一下成本。这主要涉及两个问题的计算：一个是蛋糕的质量，由此可以确定需要多少鸡蛋和面粉；另一个是蛋糕表面积(底面除外)，由此确定需要多少奶油。【实验方案】首先分析一个圆盘形的单层蛋糕，如图8-28所示。图8-28单层蛋糕绕水平中心轴旋转而成，若高为(m)，半径为r(m)，密度为(kg/m3)，则蛋糕的质量(kg)和表面积(m2)为如果蛋糕是双层圆盘的，如图8-29所示。图8-29双层蛋糕绕水平中心轴旋转而成，每层高为H/2，下层蛋糕半径为r1，上层蛋糕半径为r2，此时蛋糕的质量和表面积为以此类推，如果蛋糕是n层的如图8-30。图8-30多层蛋糕每层高为H/n，半径分别为r1，r2，……，rn，则蛋糕的质量和表面积为事实上，蛋糕边缘圆盘半径(0<h<1)那么当n→∞，H=1时此时，数学家的生日蛋糕问题就转化为求上面两个数值积分。【实验过程】symshr=2-(exp(2*h)+exp(-2*h))/5;quadl('pi*(2-(exp(2*h)+exp(-2*h))/5).^2',0,1)结果：ans=5.4171在MATLAB中输入：r0=subs(r,h,0)输出结果：r0=8/5quadl('2*pi*(2-(exp(2*h)+exp(-2*h))/5)',0,1)+pi*r0^2输出结果：ans=16.051201求得该数学家的生日大蛋糕的质量和表面积为W=5.4171(kg)，S=16.0512(m2)实验二：求侧面积为常数，体积最大的长方体体积。设长方体的长、宽、高分别为x,y,z(x>0,y>0,z>0)，体积为V，则V=f(x,y,z)=xyz。约束条件为：（x,y,z）=2(yz+zx+xy)-6=0。symsxyzlamdaa;L=x*y*z+lamda*(2*y*z+2*z*x+2*x*y-6*a^2);Lx=diff(L,'x');Ly=diff(L,'y');Lz=diff(L,'z');Llamda=diff(L,'lamda');[lamdaxyz]=solve(Lx,Ly,Lz,Llamda)运行结果：lamda=a/4-a/4x=-aay=-aaz=-aa在MATLAB中输入运行结果：V=-a^3a^3以上结果中出现的负根不在取值范围内舍去。因侧面积固定的长方体的最大体积客观存在，故当x=y=z=a时，长方体的体积最大，且最大值为。实验三：通信卫星在地面上的覆盖面积将一颗通信卫星送入太空，使该卫星轨道位于地球赤道平面内，卫星运行的角速率与地球自传的角速率（）相同时成为同步卫星。设卫星距地面的最低高度为，试计算卫星所覆盖的地球面积S。图8-31通讯卫星覆盖地面剖面图【实验方案】将地球视为球体（地球半径为），以球心为原点建立如：图8-31所示的坐标系。因上半球面方程为故被卫星覆盖的地表表面积为其中，为上半球面上被半顶角为的圆锥所截的曲面部分。所以卫星的覆盖面积为其中D:注意到于是D：利用极坐标变换，求得：当，是，S=2.1694e+008。【例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的输入输出模型中的一个基本假定是：对于每个部门，存在着一个在Rn空间单位消耗列向量vi，它表示第i个部门每产出一个单位（比如100万美金）产品，由本部门和其他各个部门消耗的百分比。在自给自足的经济中，这些列向量中所有元素的总和应该为1。把这n个vi，并列起来，它可以构成一个n×n的系数矩阵，可称为内部需求矩阵V。举一个最简单的例子，假如一个自给自足的经济体由三个部门组成，它们是煤炭业、电力业和钢铁业。它们的单位消耗列向量和销售收入列向量p如下9-1所示。表9-1单位消耗列向量和销售收入列向量p由下列部门购买每单位输出的消耗分配销售价格p（收入）煤炭业电力业钢铁业煤炭业0.00.40.6pc电力业0.60.10.2pe钢铁业0.40.50.2ps 如果电力业产出了100个单位的产品，有40个单位会被煤炭业消耗，10个单位被自己消耗，而被钢铁业消耗的是50个单位，各行业付出的费用为：这就是内部消耗的计算方法，把几个部门都算上，可以写出：各部门消耗成本==销售收入=其中： 于是总的价格平衡方程可以写成为：p–

Vp=0(I–V)p=0此等式右端常数项为零，是一个齐次方程。它有非零解的条件是系数行列式等于零，或者用行阶梯简化来求解。用MATLAB语句写出其解的表示式：>>V=[0.0.40.6;0.60.10.2;0.40.50.2];>>U0=rref([[eye(3)-V],zeros(3,1)])U0=10-31/33001-28/3300000这个结果是合理的，简化行阶梯形式只有两行，说明[I-V]的秩是2，所以它的行列式必定为零。由于现在有三个变量，只有两个方程，必定有一个变量可以作为自由变量。记住U0矩阵中各列的意义，它们分别是原方程中pc，pe，ps，的系数，所以简化行阶梯矩阵U0表示的是下列方程：这里取ps为自由变量，所以煤炭业和电力业的价格应该分别为钢铁业价格的0.94和0.85倍。如果钢铁业产品价格总计为100万元，则煤炭业的产品价格总计为94万，电力业的价格总计为85万。【例9-40】厂址选择问题。考虑A、B、C三地，每地都出产一定数量的原料也消耗一定数量的产品（见下表）。已知制成每吨产品需3吨原料，各地之间的距离为：A—B：150km，A—C：100km，B—C：200km。假定每万吨原料运输1km的运价是5000元，每万吨产品运输1km的运价是6000元。由于地区条件的差异，在不同地点设厂的生产费用也不同。问究竟在哪些地方设厂，规模多大，才能使总费用最小？另外，由于其它条件限制，在B处建厂的规模（生产的产品数量）不能超过5万吨。表9-2A、B、C三地出产原料、消耗产品情况表地点年产原料（万吨）年销产品（万吨）生产费用（万元/万吨）A207150B1613120C240100解：令为由i地运到j地的原料数量（万吨），为由i地运到j地的产品数量（万吨），i，j=1，2，3（分别对应A、B、C三地）。根据题意，可以建立问题的数学模型（其中目标函数包括原料运输费、产品运输费和生产费用（万元））：minsub.to在MATLAB中实现：>>f=[75;75;50;50;100;100;150;240;210;120;160;220];>>A=[1-11-100330000-11001-100330000-11-11000033000000001100];>>b=[20;16;24;5];>>Aeq=[000000101010000000010101];>>beq=[7;13];>>lb=zeros(12,1);>>[x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb)Optimalsolutionfound.x=010000700508fval=3485exitflag=1output=包含以下字段的struct:iterations:5constrviolation:0message:'Optimalsolutionfound.'algorithm:'dual-simplex'firstorderopt:*lambda=包含以下字段的struct:lower:[12×1double]upper:[12×1double]eqlin:[2×1double]ineqlin:[4×1double]因此，要使总费用最小，需要B地向A地运送1万吨原料，A、B、C三地的建厂规模分别为7万吨、5万吨、8万吨。最小总费用为3485万元。【例10-49】旱灾土地总面积问题。下表10-6是从1990年至2010年全国因干旱而受灾的土地总面积（单位：千公顷）数。（数据来源于全国统计局官网）表10-6各个年份全国受灾土地总面积年份受灾面积年份受灾面积年份受灾面积199018175199733516200417253199124914199814236200516028199232981199930156200620738199321097200040541200729386199430423200138472200812137199523455200222124200929259199620152200324852201013259解决以下问题：（1）计算所给样本的均值与标准差；（2）检验在显著水平为0.05的情况下，全国每年因干旱而受灾的土地总面积是否服从正态分布；（3）如果服从正态分布，用极大似然估计法对未知参数μ和σ作出估计；（4）若年受旱灾总面积大于35000千公顷即为重灾年，根据估计出的μ值和σ值，计算当年为重灾年的概率。分析问题：这是一个样本均值和标准差的计算以及正态性检验和计算的一系列问题。对于此类问题可以应用数学软件MATLAB进行处理，应用MATLAB可以很容易的计算出均值及标准差，此外，采用Jarque-Beran检验即可知道其是否服从正态分布，并估计出总体的均值μ和标准差σ。解决问题：下面计算样本的均值和标准差>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>[h,stats]=cdfplot(X)h=Line-属性:Color:[00.44700.7410]LineStyle:'-'LineWidth:0.5000Marker:'none'MarkerSize:6MarkerFaceColor:'none'XData:[1×44double]YData:[1×44double]ZData:[1×0double]显示所有属性stats=包含以下字段的struct:min:12137max:40541mean:2.4436e+04median:23455std:8.1234e+03图10-17样本经验分布函数图从输出结果可看出，样本的最小值为12137，最大值为40541，中值为23455，均值为24436，标准差为8123.4。验证下面检验其是否服从正态分布。MATLAB程序代码如下>>clearall;>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>normplot(X);>>[h,P,Jbstat,CV]=jbtest(X,0.05)>>title('正态概率图')>>xlabel('数据');>>ylable('概率')运行程序后，输出如下：h=0P=0.4574Jbstat=0.9230CV=3.8801图10-18正态概率图由输出结果h=0且Jbstat<CV可得出结论，在置信度α=0.05下，受灾面积（原始数据）服从正态分布，且在正态概率图中，各点均落在直线两侧，也可说明这一结论是成立的。再用极大似然估计法对未知参数μ和σ作出估计：MATLAB程序代码如下>>clearall;>>X=[181752491732981210973042323455201523351614236301564054138472221242485217253160282073829386121372925913259];>>phat=mle(X,'distribution','norm','alpha',0.05)运行程序后，输出如下phat=1.0e+04*2.44360.7928即受灾面积的μ估计值为24436，σ估计值为7928。最后根据估计出的μ值和σ值，计算出每年的受灾面积大于35000千公顷的概率：MATLAB程序代码如下>>clearall;p=normspec([35000inf],24436,7928)运行程序后，输出如下p=0.0913图10-19密度函数图根据输出结果可知，为重灾年的概率为0.0913.学习总结：通过对1990年至2010年全国因干旱而受灾的土地总面积的分析，我们得出这些数据服从正态分布。运用MATLAB程序，得出年均受灾土地总面积为24436公顷，根据图表可清晰地看出每年受灾总面积的分布状况，可以根据对这些数据的具体分布，采取相应的措施，从而最大程度的减小受灾。经过此次对实际问题的解决，让我们共同认识到概率统计的知识在我们的生活中无处不在，概率论在我们学习和生活中的应用也给人们带来了极大地便利。在对数据处理的过程中，对于很多数学工具的应用，如MATLAB等数学软件可让数据处理变得更加简单，从而引导我们更深层次的去探讨数学问题。在小组合作中让我们体会到小组分工合作的重要性，让我们受益匪浅。【