国赛培训：MatLab基础及在数学建模中的应用

2012年数学中国国赛培训讲座Matlab的基础及数学建模中的应用周吕文:[e-mail]中国科学院力学研究所2012年8月1第一部分MatLab基础1简单介绍MATLAB是MatrixLaboratory“矩阵实验室”的缩写。MatLab语言是由美国的CleverMoler博士于1980年开发的，初衷是为解决“线性代数”课程的矩阵运算问题。1984年由美国MathWorks公司推向市场，历经十多年的发展与竞争，现已成为国际公认的最优秀的工程应用开发环境。MATLAB功能强大、简单易学、编程效率高，深受广大科技工作者的欢迎。在数学建模竞赛中，由于只有短短的三到四天，而论文的评判不仅注重计算的结果更注重模型的创造性等很多方面，因此比赛中把大量的时间花费在编写和调试程序上只会喧宾夺主，是很不值得的。使用MATLAB可以很大程度上的方便计算、节省时间，使我们将精力更多的放在模型的完善上，所以是较为理想的。这里快速的介绍一下MATLAB与数学建模相关的基础知识，并列举一些简单的例子，很多例子都是源于国内外的数学建模赛题。希望能帮助同学们在短时间内方便、快捷的使用MATLAB解决数学建模中的问题。当然要想学好MatLab更多的依赖自主学习，一个很好的学习MatLab的方法是查看MatLab的帮助文档：如果你知道一个函数名，想了解它的用法，你可以用'help'命令得到它的帮助文档：helpfunctionname如果你了解含某个关键词的函数，你可以用'lookfor'命令得到相关的函数：2基本命令与函数基本运算变量的赋值实数赋值x=5;复数赋值x=5+10j;(或x=5+10i)向量的一般值方法行向量赋值：x=[123];(或x=[1,2,3])列向量赋值：y=[1;2;3];矩阵的赋值：x=[123;456;789];:号的用法:x=1:10/1:2:10常用矩阵（zerosoneseye）n行m列0矩阵：x=zeros(n,m);n行m列1矩阵：x=ones(n,m);n阶的单位阵：y=eye(n);rand矩阵行列操作A=[123;456;789]A=1234567892x=A(1,3)%取第一行的第三列元素A(1,3)=8;x=3y=A(2,:)%取第二行元素A(2,:)=[];y=456z=A(1:2,1:3)%取一到二行，一到三列的元素Z=123456矩阵运算（+-*/\^/‘）数组运算（点运算）A=[12;34]A=1234B=A*AB=7101522C=A.*AC=14916逻辑运算A=rand(5)A=0.15760.14190.65570.75770.70600.97060.42180.03570.74310.03180.95720.91570.84910.39220.27690.48540.79220.93400.65550.04620.80030.95950.67870.17120.0971A(A0.5)=1;A(A~=1)=0A=00111100101110001110111003常用数学函数sinsumlimitasinsortintexproundsolvelogfixdsolvelog10ceildiffsqrtfloordetabsrandsignflipud/fliplr固定变量pii,jinf–infNaNans自定义函数M文件文件名一律以字母打头函数名不能与系统函数重名function[output1,…]=functionname(input1,…)%commentofthisfunctionMatLabcommand1;MatLabcommand2;……[例2.1]做一个将角度转为孤度的函数functionRadians=deg2rad(Degrees)%DEG2RADconvertdegreestoradians%%USAGE:Radians=deg2rad(Degrees)%Degrees=[degrees,minutes,seconds]%%zhoulvwen.2011/7/25Radians=pi*Degrees*[1/180;1/180/60;1/180/60/60];4[例2.2]灰色G(1,1)模型function[q,e,E]=GM(Q,t)%E=contraryerror%e=absoluteerrorQ1=cumsum(Q);%一阶累加B1=Q1';B1(1)=[];B2=Q1';B2(end)=[];B=[-0.5*(B1+B2)ones(length(B1),1)];%构造矩阵BXn=Q';Xn(1)=[];%构造矩阵Xnc=(B'*B)^-1*B'*Xn;a=c(1);u=c(2)%计算a和uq1=(Q(1)-u/a)*exp(-a.*t)+u/a;%计算出qlq2=diff(q1);q=[q1(1)q2];%还原成qe=Q-q(1:length(Q));%算绝对误差E=e./Q;%算相对误差2流程控制语句forendifelseend/ifelseifelseend[例2.3]计算10以内的奇数和%Sumoftheoddnumbersbetween1and10Sum=0;fori=1:10ifmod(i,2)Sum=Sum+i;endendSum3二/三维作图简单作图Plotxlabelylabelxtickytickxticklabelyticklabeltextlegendaxisylimylimgridon/offgridminorsubplotlinewidthmakersize5[例3.1]在一个坐标系里绘制sin和cos曲线x=-pi:0.05:pi;y1=sin(x);y2=cos(x);plot(x,y1,'-bs',x,y2,'-m^');title('sinxandcosx');xlabel('x');ylabel('y');text(0,0,'zero');legend('sin','cos')gridonset(gca,'xtick',[-pi:pi/2:pi])xlim([-pi,pi])set(gca,’xticklabel’,)set(gca,'xticklabel',{'-pi','-pi/2','0','pi/2','pi'})set(gca,'XTickLabel',{'-p','-p/2','0','p/2','p'},'FontName','Symbol')-4-3-2-101234-1-0.500.51sinxandcosxxyzerosincos作图函数fplotfplot('sin(x)',[-pi,pi])polart=0:.01:2*pi;polar(t,sin(2*t).*cos(2*t),'-r')barbar(1:4,[35,23,9,20])piecontour[x,y,z]=peaks;contour(x,y,z)quiverquiver(x,y,px,py)imageplot3meshgrid[x,y]=meshgrid(1:5,1:5)meshsurfcontour3contour3(peaks)meshcsurfc[x,y]=meshgrid(1:5,1:5)6[例3.2]用plot3作三维螺旋线z=0:pi/50:10*pi;x=sin(z);y=cos(z);plot3(x,y,z)title('Helix')xlabel('sint(t)')ylabel('cos(t)')zlabel('t')text(0,0,0,'Origin')gridon-101-101010203040sint(t)OriginHelixcos(t)t[例3.3]用mesh作peaks函数的网格图。[X,Y,Z]=peaks(30);mesh(X,Y,Z)gridonxlabel('x-axis')ylabel('y-axis')zlabel('z-axis')title('MESHofPEAKS')-3-2-10123-3-2-10123-8-6-4-20246810x-axisMESHofPEAKSy-axisz-axis7第二部分例题精讲1人口问题表1是1790年至1900年美国人口数(单位:百万)，请用指数增长模型拟合美国人口变化。表1.1790年至1900年美国人口数（单位：百万）17903.9184017.1189062.918005.3185023.2190076.018107.2186031.418209.6187038.6183012.9188050.2参考程序：prediction_of_population.m%translatethequstionintopolynomialfit,thenwecanusethe%polyfitfuction%P=x0*exp(r*t)=lnP=lnx0+r*t%\|/\|/\|/\|/%Y=a2+a1*xt=1790:10:1900;population=[3.95.37.29.612.917.123.231.438.650.262.976.0];Y=log(population);X=t;a=polyfit(X,Y,1);x0=exp(a(2))r=a(1)ti=1790:1950;poi=x0*exp(r*ti);plot(t,population,'k.',ti,poi,'b','markersize',18,'linewidth',2);legend('Actualdata','Fitresult')xlabel('Year');ylabel('Population')82图像处理在Matlab中,图可以表示成矩阵,矩阵也可能表示成图。对图像的处理就可转化为对矩阵的处理（见附件1）。设计算法求出任一矩阵的边界矩阵（算法应具有一般性）。边界矩阵是要找出图形的边界。下面为边界矩阵的定义：BB的边界矩阵图1边界矩阵的定义参考程序：Boundary.m,main_of_boundary.mfunctionboundary=Boundary(M)[I,J]=size(M);i=2:I-1;j=2:J-1;SUM=M(i-1,j)+M(i+1,j)+M(i,j-1)+M(i,j+1);%Sumoffourneighboursm=M(i,j);m(SUM==4)=0;%ifallof(i,j)'sfourneighbours==1thenm(i,j)=0M(i,j)=m;boundary=M;M=load('A.txt');boundary=Boundary(M);subplot(1,2,1)image(M*255)%showtheMmatrixsubplot(1,2,2)image(boundary*255)%showtheboundarymatrix0010001110111111111001100001000101010001100100110093矩阵构造用元胞自动机模型研究收费站设置多少收费亭最佳时,需要根据给定的高速公路(车道数L),和收费亭数目B以及与收费站结构有关的参数生成相应的格子(矩阵).每个格子可能的状态有三种:用1表示车辆,0表示空位,-888表示不可进入区域.图2收费站的离散化参考程序：create_plaza.mfunction[plaza,v,time]=create_plaza(B,L,plazalength)plaza=zeros(plazalength,B+2);%1=car,0=empty,-888=forbidv=zeros(plazalength,B+2);%velocityofautomata(i,j),ifitexiststime=zeros(plazalength,B+2);%costtimeofautomata(i,j)ifitexistsplaz