(2)MATLAB应用实例分析

1Matlab应用例题选讲仅举一些运用MATLAB的例子，这些问题在数学建模中时常遇到，希望能帮助同学们在短时间内方便、快捷的使用MATLAB解决数学建模中的问题，并善用这一工具。常用控制命令：clc：%清屏；clear：%清变量；save：%保存变量；load：%导入变量一、利用公式直接进行赋值计算本金P以每年n次，每次i%的增值率（n与i的乘积为每年增值额的百分比）增加，当增加到r×P时所花费的时间T为：（利用复利计息公式可得到下式）)01.01ln(ln)01.01(inrTiPPrnT（12,5.0,2nir）MATLAB的表达形式及结果如下：r=2;i=0.5;n=12;%变量赋值T=log(r)/(n*log(1+0.01*i))计算结果显示为：T=11.5813即所花费的时间为T=11.5813年。分析：上面的问题是一个利用公式直接进行赋值计算问题，实际中若变量在某个范围变化取很多值时，使用MATLAB，将倍感方便，轻松得到结果，其绘图功能还能将结果轻松的显示出来，变量之间的变化规律将一目了然。若r在[1,9]变化，i在[0.5,3.5]变化；我们将MATLAB的表达式作如下改动，结果如图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'))1234567890510152025303540rT0.511.522.533.5图12从图1中既可以看到T随r的变化规律，而且还能看到i的不同取值对T—r曲线的影响（图中的六条曲线分别代表i的不同取值）。二、已知多项式求根已知多项式为9620011610311023456xxxxxxh，求其根。分析：对多项式求根问题，我们常用roots()函数。MATLAB的表达形式及结果如下：h=roots([1-1031-10-116200-96])%中括号内为多项式系数由高阶到常数。计算结果显示为（其中i为虚数单位）：h=-2.00004.00003.00002.0000+0.0000i2.0000-0.0000i1.0000如果已知多项式的根，求多项式，用poly()函数。对上面得到的h的值求多项式，其MATLAB的表达形式及结果如下：h=[-2.00004.00003.00002.0000+0.0000i2.0000-0.0000i1.0000];c=poly(h)计算结果显示为：c=1-1031-10-116200-96三、方程组的求解求解下面的方程组：1229447535.768321321321xxxxxxxxx分析：对于线性方程组求解，常用线性代数的方法，把方程组转化为矩阵进行计算。baxbax1bax\MATLAB的表达形式及结果如下：a=[816;357;492];%建立系数矩阵b=[7.5;4;12];%建立常数项矩阵x=a\b%求方程组的解计算结果显示为：x=1.29310.8972-0.6236四、数据拟合与二维绘图在数学建模竞赛中，我们常会遇到这种数据表格问题，如果我们仅凭眼睛观察，很难看到其中的规律，也就更难写出有效的数学表达式从而建立数学模型。因此可以利用MATLAB的拟合函数，即polyfit()函数，并结合MATLAB的绘图功能（利用plot()函数），得到直观的表示。3例：在化学反应中，为研究某化合物的浓度随时间的变化规律，测得一组数据如下表：分析：MATLAB的表达形式如下：t=[1:16];%数据输入y=[46.488.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];plot(t,y,'o')%画散点图p=polyfit(t,y,2)%二次多项式拟合holdonxi=linspace(0,16,160);%在[0,16]等间距取160个点yi=polyval(p,xi);%由拟合得到的多项式及xi，确定yiplot(xi,yi)%画拟合曲线图执行程序得到图2；02468101214164567891011图2显示的结果为p=-0.04451.07114.3252p的值表示二阶拟合得到的多项式为：y=-0.0445t2+1.0711t+4.3252下面是用lsqcurvefit()函数，即最小二乘拟合方法的Matlab表达：t=[1:16];y=[46.488.49.289.59.79.861010.210.3210.4210.510.5510.5810.6];x0=[0.1,0.1,0.1];zuixiao=inline('x(1)*t.^2+x(2)*t+x(3)','x','t');x=lsqcurvefit(zuixiao,x0,t,y)%利用最小二乘拟合其显示的结果为：x=-0.04451.07114.3252可以看出其得到的结果与polyfit函数的结果相同。这说明在多项式拟合问题上这两个函数的效果是相同的。T（分）12345678y46.48.08.49.289.59.79.86T（分）910111213141516y1010.210.3210.4210.510.5510.5810.64下面的一个例子将体现lsqcurvefit()函数的优势。例2:在物理学中，为研究某种材料应力与应变的关系，测得一组数据如下表：应力σ925112516252125262531253625应变ε0.110.160.350.480.610.710.85如果假定应力与应变有如下关系（σ为应力值，ε为应变值）：ε=a+blnσ试计算a、b的值。MATLAB的表达形式如下：x=[925,1125,1625,2125,2625,3125,3625];y=[0.11,0.16,0.35,0.48,0.61,0.71,0.85];plot(x,y,'o')[p,resid1]=polyfit(x,y,2)holdonxi=linspace(700,3700,3000);yi=polyval(p,xi);plot(xi,yi)x0=[0.1,0.1];fff=inline('a(1)+a(2)*log(x)','a','x');[a,resid2]=lsqcurvefit(fff,x0,x,y)plot(xi,fff(a,xi),'r')执行程序得到图3，图中蓝色曲线为利用polyfit()函数得到的曲线，红色曲线为利用lsqcurvefit()函数得到的曲线；5001000150020002500300035004000-0.100.10.20.30.40.50.60.70.80.9其显示的结果为：p=-0.00000.0004-0.2266resid1=R:[3x3double]df:4normr:0.03315a=-3.58100.5344resid2=0.0064其中a的值代表利用lsqcurvefit()函数得到的关系为：ε=-3.5810+0.5344σresid1、resid2分别代表运用polyfit()函数、lsqcurvefit()函数得到的残差。可以看出利用lsqcurvefit()函数残差更小，即得到了更好的拟合效果。在数学建模的实际问题中，如果问题的机理不明，我们只能采用polyfit()函数，即多项式拟合的方法，以获得近似的数据描述函数；但如果通过分析，可以得到一些机理，那么采用最小二乘的方法将得到更好的效果，而且得到的拟合函数也更有意义。五、隐函数的图形绘制plot()只能绘制显函数图形，对于形如0)sin()1ln(ln1xxyyy的复杂隐函数，很难转化为显函数并利用plot()函数绘制图形，这时就可以用ezplot()函数直接绘制其曲线。MATLAB的表达形式如下：ezplot('1/y-log(y)+log(-1+y)+x-sin(x)')执行程序得到图5-6-4-20246-6-4-20246xy1/y-log(y)+log(-1+y)+x-sin(x)=0-0.8-0.6-0.4-0.200.20.40.60.8-0.8-0.6-0.4-0.200.20.4xyx=sin(3t)cos(t),y=sin(3t)sin(t)图5图6如果是形如下面的参数方程),0(,sin3sin,cos3sintttyttx,同样可以利用ezplot()函数绘制其曲线。MATLAB的表达形式如下：ezplot('sin(3*t)*cos(t)','sin(3*t)*sin(t)',[0,pi])执行程序得到图6。六、三维图形绘制假设有一个时间向量t，对该向量进行下列运算则可以构成三个坐标值向量tztytx,cos,sin对于上面的方程可以利用ezplot3()函数或plot3()函数绘制三维曲线。这里仅列举ezplot3()函数的使用。MATLAB的表达形式如下：ezplot3('sin(t)','cos(t)','t',[0,6*pi])执行程序得到图7：绘制下述曲面：20,10),3cos(),(3rrrz其中MATLAB的表达形式如下：6nr=12;nth=50;r=linspace(0,1,nr);theta=linspace(0,2*pi,nth);[R,T]=meshgrid(r,theta)x=cos(theta')*r;y=sin(theta')*r;surf(x,y,R.^3.*cos(3*T))执行程序得到图8。-1-0.500.51-1-0.500.5105101520xx=sin(t),y=cos(t),z=tyz图7图8除了surf()函数还有surfc()、surfl()、mesh()、waterfall()函数也用于曲面的绘制，具体效果如图9所示，可以针对自己的需要选取适合的曲面绘制函数。图9MATLAB作图之三维绘图示例山体绘制%三维绘图示例山体绘制%mesh函数演示x=1.0:0.1:2.0;y=2.0:0.1:3.0;[X,Y]=meshgrid(x,y);z=[5.115.135.145.135.095.044.984.934.894.854.855.395.495.515.465.325.144.944.744.594.494.485.615.775.815.715.515.234.904.594.364.214.195.735.925.975.865.625.274.884.514.234.054.035.745.925.975.867.625.274.884.514.214.044.025.635.795.846.7410.539.238.914.594.334.184.165.425.535.565.497.355.164.934.734.554.454.445.145.185.195.1711.125.054.974.904.844.814.8074.484.804.794.824.874.945.025.105.165.195.204.564.454.434.494.644.845.065.285.455.555.564.364.194.164.254.474.765.095.415.665.815.83];figure(1)%画网格图mesh(X,Y,z);colormap([010]);xlabel('x轴');ylabel('y轴');zlabel('z轴');%画表面图figure(2)surf(X,Y,z);colormap([100]);xlabel('x轴');ylabel('y轴');zlabel('z轴');11.2