2011年数学中国国赛培训讲座 Matlab的基础及数学建模中的应用周吕文:zhou.lv.wen@大连大学数学建模工作室&中国科学院力学研究所2011年7月第一部分 MatLab基础1 简单介绍MATLAB是Matrix Laboratory“矩阵实验室”的缩写。

MatLab语言是由美国的Clever Moler博士于1980年开发的，初衷是为解决“线性代数”课程的矩阵运算问题。

1984年由美国 MathWorks公司推向市场，历经十多年的发展与竞争，现已成为国际公认的最优秀的工程应用开发环境。

MATLAB功能强大、简单易学、编程效率高，深受广大科技工作者的欢迎。

在数学建模竞赛中，由于只有短短的三到四天，而论文的评判不仅注重计算的结果更注重模型的创造性等很多方面，因此比赛中把大量的时间花费在编写和调试程序上只会喧宾夺主，是很不值得的。

使用MATLAB 可以很大程度上的方便计算、节省时间，使我们将精力更多的放在模型的完善上，所以是较为理想的。

这里快速的介绍一下MATLAB与数学建模相关的基础知识，并列举一些简单的例子，很多例子都是源于国内外的数学建模赛题。

希望能帮助同学们在短时间内方便、快捷的使用MATLAB 解决数学建模中的问题。

当然要想学好MatLab更多的依赖自主学习，一个很好的学习MatLab的方法是查看MatLab的帮助文档：z如果你知道一个函数名，想了解它的用法，你可以用'help'命令得到它的帮助文档：>>help functionnamez如果你了解含某个关键词的函数，你可以用'lookfor'命令得到相关的函数：2 基本命令与函数基本运算z变量的赋值实数赋值>> x=5;复数赋值>> x=5+10j; (或>>x=5+10i)z向量的一般值方法行向量赋值：>>x=[1 2 3]; (或x=[1, 2 ,3])列向量赋值：>>y=[1;2;3];矩阵的赋值：>>x=[1 2 3; 4 5 6; 7 8 9];z常用矩阵（zeros ones eye）n行m列0矩阵：>>x=zeros(n,m);n行m列1矩阵：>>x=ones(n,m);n 阶的单位阵：>>y=eye(n);z矩阵行列操作>> A=[1 2 3;4 5 6;7 8 9]A=1 2 34 5 67 8 9>>x=A(1,3) %取第一行的第三列元素x=3>>y=A(2, :) %取第二行元素y=4 5 6>>z=A(1:2,1:3)%取一到二行，一到三列的元素Z=1 2 34 5 6z矩阵运算（+ - * / \ ^ / ‘）z数组运算（点运算）>>A=[1 2; 3 4]A=1 23 4>>B=A*AB=7 1015 22>>C=A.*AC=1 49 16常用数学函数sin sum limit asin sort intsolve exp roundlog fix dsolve log10 ceil diff sqrt floor det abs randsign flipud/fliplr固定变量pii,jinf –infNaNans自定义函数z M文件文件名一律以字母打头z函数名不能与系统函数重名function [output 1, …] = functionname(input1, …)% comment of this functionMatLab command 1;MatLab command 2;… …[例2.1]做一个将角度转为孤度的函数function% DEG2RAD convert degrees to radians%% USAGE: Radians = deg2rad(Degrees)% Degrees =[degrees, minutes, seconds]%% zhou lvwen. 2011/7/25Radians = pi * Degrees * [1/180; 1/180/60; 1/180/60/60];[例2.2]灰色G(1,1)模型function [q,e,E] = GM(Q,t)% E = contrary error% e = absolute errorQ1 = 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 流程控制语句z for endz if else end / if elseif else end[例2.3]计算10以内的奇数和% Sum of the odd numbers between 1 and 10Sum = 0;for i = 1:10if mod(i,2)Sum =Sum + i;endendSum3 二/三维作图简单作图z Plotz Xlabel ylabel xtick ytick xticklabel yticklabel z text legendz axis ylim ylimz grid on/off grid minorz subplotz linewidth makersize作图函数fplot polar barpie contour quiver image plot3meshgridmesh surf contour3meshc surfc[例3.2]用plot3作三维螺旋线x = -pi:0.05:pi;y1 = sin(x);y2 = cos(x);plot(x,y1,'-bs',x,y2,'-m^');title('Helix');xlabel('x');ylabel('y');text(0,0,'zero');legend('sin','cos')grid on第二部分例题精讲1人口问题表1是1790年至1900年美国人口数(单位:百万)，请用指数增长模型拟合美国人口变化。

表1. 1790 年至1900年美国人口数（单位：百万）1790 3.9 1840 17.1 1890 62.91800 5.3 1850 23.2 1900 76.01810 7.2 1860 31.41820 9.6 1870 38.61830 12.9 1880 50.2参考程序：prediction_of_population.m% translate the qustion into polynomial fit,then we can use the% polyfit fuction% P = x0 * exp(r*t) => lnP = lnx0 + r * t% \|/ \|/ \|/ \|/% Y = a2 + a1* xt=1790:10:1900;population=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.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('Actual data', 'Fit result')xlabel('Year');ylabel('Population')2图像处理在Matlab 中,图可以表示成矩阵,矩阵也可能表示成图。

对图像的处理就可转化为对矩阵的处理（见附件1）。

设计算法求出任一矩阵的边界矩阵（算法应具有一般性）。

边界矩阵是要找出图形的边界。

下面为边界矩阵的定义：B B 的边界矩阵图1边界矩阵的定义参考程序：Boundary.m, main_of_boundary.m function boundary=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);%Sum of four neighbours m=M(i,j);m(SUM==4)=0; % if all of (i,j)'s four neighbours==1 then m(i,j)=0 M(i,j)=m; boundary=M;'A.txt'boundary=Boundary(M); subplot(1,2,1)image(M*255) % show the M matrix subplot(1,2,2)image(boundary*255) % show the boundary matrix00 1 0 00 1 1 1 01 1 1 1 11 1 1 1 00 1 1 0 000100010101000110010011003矩阵构造用元胞自动机模型研究收费站设置多少收费亭最佳时,需要根据给定的高速公路(车道数L),和收费亭数目B以及与收费站结构有关的参数生成相应的格子(矩阵).每个格子可能的状态有三种:用1表示车辆,0表示空位,-888表示不可进入区域.图2 收费站的离散化参考程序：create_plaza.mfunctionplaza = zeros(plazalength,B+2); % 1 = car, 0 = empty, -888 = forbidv = zeros(plazalength,B+2); % velocity of automata (i,j), if it existstime = zeros(plazalength,B+2); % cost time of automata (i,j) if it existsplaza(1:plazalength,[1,2+B]) = -888;%left: angle of width decline for boundariestoptheta = 1.3;bottomtheta = 1.2;for col = 2:ceil(B/2-L/2) + 1for row = 1:(plazalength-1)/2 - floor(tan(toptheta) * (col-1))plaza(row, col) =-888;endfor row = 1:(plazalength-1)/2 - floor(tan(bottomtheta) * (col-1))plaza(plazalength+1-row, col) =-888;endendfac = ceil(B/2-L/2)/floor(B/2-L/2);%right: angle of width decline for boundariestoptheta = atan(fac*tan(toptheta));bottomtheta = atan(fac*tan(bottomtheta));forfor row = 1:(plazalength-1)/2 - floor(tan(toptheta) * (col-1))plaza(row,B+3-col) =-888;endfor row = 1:(plazalength-1)/2 - floor(tan(bottomtheta) * (col-1))plaza(plazalength+1-row,B+3-col) = -888;endend%% Alternative method %%% topgap = 5;% bottomgap = 1;% plaza = zeros(plazalength,B+2);% plaza(1:plazalength,[1,2+B]) = -888;% if mod(B-L,2)==0% for col = 2:B/2 - L/2 + 1% for row = 1:(plazalength-1)/2 - topgap * (col-1)% plaza(row,[col, B+3-col]) = -888;% end% for row =(plazalength+3)/2 + bottomgap*(col-1):plazalength % plaza(row,[col, B+3-col]) = -888;% end% end% else% plaza(1:plazalength, B+3) =-888;% for col = 2:(B+1)/2 - L/2 + 1% for row = 1:(plazalength-1)/2 - topgap * (col-1)% plaza(row, [col, B+4-col]) = -888;% end% for row =(plazalength+3)/2 + bottomgap*(col-1):plazalength % plaza(row, [col, B+4-col]) = -888;% end% end% end%%4地面搜索地震后，有一块km 5050×的地面须要搜索，以抢救伤员。