《数学建模软件求解入门.ppt》由会员分享，可在线阅读，更多相关《数学建模软件求解入门.ppt（83页珍藏版）》请在三一办公上搜索。
1、Matlab与数学建模,Matlab简介,1997年仲春，MATLAB5.0版问世，紧接着是5.1、5.2，以及和1999年春的5.3版。与 4.x相比，现今的MATLAB拥有更丰富的数据类型和结构、更友善的面向对象、更加快速精良的图形可视、更广博的数学和数据分析资源、更多的应用开发工具。发展至今，Matlab7.0版已经问世,Matlab简介,Matlab简介,Matlab软件主要由三部分组成：Matlab主包，Simulink和工具箱,Matlab软件的组成,第1章 矩阵及其基本运算,1实数值矩阵输入MATLAB的强大功能之一体现在能直接处理向量或矩阵。当然首要任务是输入待处理的向量或矩阵
2、。不管是任何矩阵（向量），我们可以直接按行方式输入每个元素：同一行中的元素用逗号（，）或者用空格符来分隔，且空格个数不限；不同的行用分号（；）分隔。所有元素处于一方括号（）内；当矩阵是多维（三维以上），且方括号内的元素是维数较低的矩阵时，会有多重的方括号。如：,1)Time=11 12 1 2 3 4 5 6 7 8 9 10Time=11 12 1 2 3 4 5 6 7 8 9 10,2)A=1 2 3；2 3 4；3 4 5A=1 2 32 3 43 4 5,3)Null_M=%生成一个空矩阵,4)A=A;1 2 3,1;2;3;4,特殊矩阵的生成,命令 全零阵函数 zeros格式 B=
3、zeros(n)%生成nn全零阵B=zeros(m,n)%生成mn全零阵B=zeros(m n)%生成mn全零阵B=zeros(size(A)%生成与矩阵A相同大小的全零阵,命令 单位阵函数 eye格式 Y=eye(n)%生成nn单位阵Y=eye(m,n)%生成mn单位阵Y=eye(size(A)%生成与矩阵A相同大小的单位阵,命令 均匀分布随机矩阵函数 rand 格式 Y=rand(n)%生成nn随机矩阵，其元素在（0，1）内Y=rand(m,n)%生成mn随机矩阵,命令 全1阵函数 ones格式 Y=ones(n)%生成nn全1阵Y=ones(m,n)%生成mn全1阵Y=ones(m n)
4、%生成mn全1阵Y=ones(size(A)%生成与矩阵A相同大小的全1阵,思考：如何产生20至50间的随机数？,向量生成方法,2.定数线性采样法,该采样法在给定的范围内确定等距离的一个样本数。该方法很实用，特别适合定义大的数组,例如：在 之间取5个点,linspace(-pi,pi,5),3.利用随机数发生器定义一维数组,例如：产生5个（0，1）均匀分布随机数，并定义一维数组x。,x=rand(1,5),Matlab常用命令,Matlab常用命令,若要检视现存於工作空间（Workspace）的变数，可键入who：who Your variables are:testfile x 这些是由使用
5、者定义的变数。若要知道这些变数的详细资料，可键入：whos Name Size Bytes Class A 2x4 64 double array B 4x2 64 double array ans 1x1 8 double array x 1x1 8 double array y 1x1 8 double array z 1x1 8 double array Grand total is 20 elements using 160 bytes 使用clear可以删除工作空间的变数：clear A A?Undefined function or variable A.,Matlab常用命令应用
6、,变量与常量变量是任何程序设计语言的基本要素之一，MATLAB语言当然也不例外。与常规的程序设计语言不同的MATLAB并不要求事先对所使用的变量进行声明，也不需要指定变量类型，MATLAB语言会自动依据所赋予变量的值或对变量所进行的操作来识别变量的类型。在赋值过程中如果赋值变量已存在时，MATLAB语言将使用新值代替旧值，并以新值类型代替旧值类型。在MATLAB语言中变量的命名应遵循如下规则：（1）变量名区分大小写。（2）变量名长度不超31位，第31个字符之后的字符将被MATLAB语言所忽略。（3）变量名以字母开头，可以是字母、数字、下划线组成，但不能使用标点。与其他的程序设计语言相同，在MA
7、TLAB语言中也存在变量作用域的问题。在未加特殊说明的情况下，MATLAB语言将所识别的一切变量视为局部变量，即仅在其使用的M文件内有效。若要将变量定义为全局变量，则应当对变量进行说明，即在该变量前加关键字global。一般来说全局变量均用大写的英文字符表示。MATLAB语言本身也具有一些预定义的变量，这些特殊的变量称为常量。,另外MATLAB有些永久常数（Permanent constants），虽然在工作空间中看不 到，但使用者可直接取用，例如：pi ans=3.1416=下表即为MATLAB常用到的永久常数。小整理：MATLAB的永久常数 i或j：基本虚数单位（即）eps：系统的浮点（F
8、loating-point）精确度 inf：无限大，例如1/0 nan或NaN：非数值（Not a number），例如0/0 pi：圆周率 p（=3.1415926.）realmax：系统所能表示的最大数值 realmin：系统所能表示的最小数值 nargin:函数的输入引数个数 nargout:函数的输出引数个数,特别注意,实现Matlab与word的无缝连接,步骤：1）在命令窗口输入notebook-setup 2)选择word版本 3）在命令窗口输入notebook 4)word中输入文字以及Matlab语句，全部选中，选中下拉菜单，点选define input cell，再选择eva
9、luate cell,作出泊松分布的密度图形x=0:20;y=poisspdf(x,5);stem(x,y),将硬盘的数据读入变量,现在有一组数据保存在three.xls文件中，要在Matlab中读入这组数据,a=xlsread(e:three.xls),如果保存在文件soil.txt中只要输入a=load(e:soil.txt),矩阵的简单运算,结果显示：A+B=9 2 74 7 105 12 8AB=-7 0-5-2-3-4-3-6 4,运算符：“”和“”分别为加、减运算符。运算规则：对应元素相加、减，即按线性代数中矩阵的“十”，“一”运算进行。,加、减运算,例A=1,1,1;1,2,3;
10、1,3,6B=8,1,6;3,5,7;4,9,2AB=A+BA-=A-B,运算符：*运算规则：按线性代数中矩阵乘法运算进行，即放在前面的矩阵的各行元素，分别与放在后面的矩阵的各列元素对应相乘并相加。1两个矩阵相乘,乘法,例X=2 3 4 5;1 2 2 1；Y=0 1 1;1 1 0;0 0 1;1 0 0；Z=X*Y,结果显示为：Z=8 5 6 3 3 3,上例中：a=2*X则显示：a=4 6 8 102 4 4 2,矩阵的数乘：数乘矩阵,维数相同的两个向量的点乘。数组乘法：A.*B表示A与B对应元素相乘。,向量的点乘（内积）,需要注意的是点乘与矩阵的乘法是有区别的,例：a=1 2 3;4
11、2 6;7 4 9b=4;1;2;x=ab则显示：x=-1.5000 2.00000.5000如果a为非奇异矩阵，则ab和b/a可通过a的逆矩阵与b阵得到：ab=inv(a)*b b/a=b*inv(a)数组除法：A./B表示A中元素与B中元素对应相除。,矩阵除法,Matlab提供了两种除法运算：左除（）和右除（/）。一般情况下，x=ab是方程a*x=b的解，而x=b/a是方程x*a=b的解。,同样此时也要注意除法和点除之间的区别,运算符：运算规则：（1）当A为方阵，P为大于0的整数时，AP表示A的P次方，即A自乘P次；P为小于0的整数时，AP表示A-1的P次方。（4）标量的数组乘方P.A，标
12、量的数组乘方定义为数组乘方：A.P：表示A的每个元素的P次乘方。,矩阵乘方,矩阵的乘方要求矩阵A是方阵,例 A=1 2 3;4 5 6;7 8 9A=1 2 3 4 5 6 7 8 9 D=det(A)D=0,运算符：运算规则：若矩阵A的元素为实数，则与线性代数中矩阵的转置相同。若A为复数矩阵，则A转置后的元素由A对应元素的共轭复数构成。若仅希望转置，则用如下命令：A.。,矩阵转置,方阵的行列式,函数 det格式 d=det(X)%返回方阵X的行列式的值,函数 inv格式 Y=inv(X)%求方阵X的逆矩阵。若X为奇异阵或近似奇异阵，将给出警告信息。例1-43 求的逆矩阵方法一A=1 2 3;
13、2 2 1;3 4 3;Y=inv(A)或Y=A(-1)则结果显示为 Y=1.0000 3.0000-2.0000-1.5000-3.0000 2.5000 1.0000 1.0000-1.0000方法二：由增广矩阵进行初等行变换B=1,2,3,1,0,0;2,2,1,0,1,0;3,4,3,0,0,1；C=rref(B)%化行最简形X=C(:,4:6)%取矩阵C中的A(-1)部分,矩阵求逆,命令 逆,子矩阵的提取和运算,使用命令：BA（v1，v2），其中，v1向量表示子矩阵要包含的行号构成的向量，v2表示要包含的列号构成的向量。若v1为：，则表示要提取所有的行，同样使用列号。关键词end表示
14、最后一行（或列）,例如：B1A(1:2:end,:),表示提取A矩阵的全部奇数行、所有列,B2=A(3,2,1,2,3,4),表示提取A矩阵3，2，1行、2，3，4列构成子矩阵。,B3=A(:,end:-1:1),表示将A矩阵左右翻转，即最后一列排在最前面。,A(:,2)=%删除第二列（：代表所有行）,例 A=1 2 3 4 5 6 7 8 9 10 11 12,我们可以对矩阵进行各种处理：,A(2,3)=5%改变位於第二列，第三行的元素值 A=1 2 3 4 5 6 5 8 9 10 11 12,B=A(2,1:3)%取出部份矩阵B B=5 6 5,A=A B%将B转置後以行向量并入A A=
15、1 2 3 4 5 5 6 5 8 6 9 10 11 12 5,A=A;4 3 2 1%加入第四行 A=1 3 4 5 5 5 8 6 9 11 12 5 4 3 2 1,A(1 4,:)=%删除第一和第四行（：代表所有列）A=5 5 8 6 9 11 12 5,数组操作函数,Repmat（D，1，3）命令：在水平方向“铺放”三个矩阵,D1 0 00 1 00 0 1,则Repmat（D，1，3）结果为,ans=1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1,Diag(A,1)命令：对矩阵A取第一上对角元素,可运行diag(A,-1
16、),查看运行结果,Triu(A)命令：取A的上三角阵；tril(A)命令：取A的下三角阵,矩阵的逻辑运算,矩阵的比较运算：用表示等于关系，用=表示大于等于关系，用=表示不等于关系等。,A=1 2 3;4,5,6;7,8,0Find(A=5),找出矩阵中元素大于等于5的下标,ans=3 5 6 8,该函数先将A按列构成列向量，然后判断哪些元素大于等于5，返回其下标。,i,j=find(A=5);i,j,同时返回行和列坐标,all(A=5),判断A矩阵的某列元素是否全部大于等于5，是返回1，否则为0,any(A=5),你能猜测表示什么含义吗？,abs(x)：纯量的绝对值或向量的长度 angle(z
17、)：复数z的相角(Phase angle)sqrt(x)：开平方 real(z)：复数z的实部 imag(z)：复数z的虚部 conj(z)：复数z的共轭复数 round(x)：四舍五入至最近整数 fix(x)：无论正负，舍去小数至最近整数 floor(x)：地板函数，即舍去正小数至最近整数 ceil(x)：天花板函数，即加入正小数至最近整数 rat(x)：将实数x化为分数表示 rats(x)：将实数x化为多项分数展开 sign(x)：符号函数(Signum function)。当x0时，sign(x)=1。rem(x,y)：求x除以y的馀数 gcd(x,y)：整数x和y的最大公因数 lcm(
18、x,y)：整数x和y的最小公倍数,MATLAB常用的基本数学函数,第2章 绘图与图形处理,命令1 plot功能 线性二维图。在线条多于一条时，若用户没有指定使用颜色，则plot循环使用由当前坐标轴颜色顺序属性（current axes ColorOrder property）定义的颜色，以区别不同的线条。在用完上述属性值后，plot又循环使用由坐标轴线型顺序属性（axes LineStyleOrder property）定义的线型，以区别不同的线条。用法 plot(X,Y)当X,Y均为实数向量，且为同维向量（可以不是同型向量），X=x(i)，Y=y(i)，则plot(X,Y)先描出点(x(i)
19、，y(i)，然后用直线依次相连。,二维图形的绘制,1 散点图的绘制,画散点图一般使用命令scatter，其调用格式为：scatter(x,y),2 二维曲线的绘制,plot(X1,Y1,LineSpec1,X2,Y2,LineSpec2)将按顺序分别画出由三参数定义Xi,Yi,LineSpeci的线条。其中参数LineSpeci指明了线条的类型，标记符号，和画线用的颜色。,例7-1 t=0:pi/20:2*pi;plot(t,t.*cos(t),-.r*)hold on plot(exp(t/100).*sin(t-pi/2),-mo)plot(sin(t-pi),:bs)hold off,绘
20、图中的一些线型和颜色,特殊标记,theta=0:0.01:6*pi;rho=5*sin(4*theta/3);polar(theta,rho),例子：试用极坐标绘制函数polar()绘制出 的极坐标曲线,分块显示的例子,关于条形图的例子,建立三维饼图,三维图形的绘制,1 三维曲线的绘制,三维曲线可以使用plot()命令，其调用格式为,plot3(x,y,z),theta=0:0.01:2*pi;x=sin(theta);y=cos(theta);z=theta;plot3(x,y,z),例如：以下的程序是在绘制什么图形？,2 三维曲面绘制,绘制三维曲面，先调用meshgrid()函数生成网格矩
21、阵数据x和y，再使用mesh()或surf()等函数进行三维图形绘制。其调用格式为,x,y=meshgrid(v1,v2);z=x.*y;mesh(x,y,z),练习：求x2+y2=z的空间曲面,要求：图像的横，纵坐标的范围在-4到4之间,1、在图形上加格栅、图例和标注,1）GRID ON:加格栅在当前图上GRID OFF:删除格栅,处理图形,2）hh=xlabel(string):在当前图形的x轴上加图例string,hh=ylabel(string):在当前图形的y轴上加图例string,hh=title(string):在当前图形的顶端上加图例string,hh=zlabel(strin
22、g):在当前图形的z轴上加图例string,例 在区间0,2*pi画sin(x)的图形，并加注图例“自变量 X”、“函数Y”、“示意图”,并加格栅.,解 x=linspace(0,2*pi,30);y=sin(x);plot(x,y)xlabel(自变量X)ylabel(函数Y)title(示意图)grid on,Matlab liti2,3）hh=gtext(string),命令gtext(string)用鼠标放置标注在现有的图上.运行命令gtext(string)时，屏幕上出现当前图形，在图形上出现一个交叉的十字，该十字随鼠标的移动移动，当按下鼠标左键时，该标注string放在当前十交叉的
23、位置.,例 在区间0,2*pi画sin(x)，并分别标注“sin(x)”cos(x)”.,解 x=linspace(0,2*pi,30);y=sin(x);z=cos(x);plot(x,y,x,z)gtext(sin(x);gtext(cos(x),Matlab liti3,返回,第三章 基本编程,循环语句,这个例子想做什么？,从1到多少的自然数和大于或等于100,循环语句,当有一个等效的数组方法来解给定的问题时，应避免用For循环。例如 for n=1:10 x(n)=sin(n*pi/10);end xx=Columns 1 through 7 0.3090 0.5878 0.8090
24、0.9511 1.0000 0.9511 0.8090 Columns 8 through 10 0.5878 0.3090 0.0000 这个例子可重写为：n=1:10;x=sin(n*pi/10)x=Columns 1 through 7 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 Columns 8 through 10 0.5878 0.3090 0.0000两种方法得出同样的结果，而后者执行更快，更直观，要求较少的输入。,循环语句,条件语句,条件语句,条件语句,例子,if x=1 y=10;elseif x-1&x1 y=0;el
25、se y=-10;end,编程如下,Matlab ifprom,条件语句,条件语句,t=-pi:0.1:pi;trigname=input(input trig function name:);switch trigname case sin plot(t,sin(t)case cos plot(t,cos(t)case tan plot(t,tan(t)otherwise breakend,Matlab switchprom,【*例8.1-1】通过M脚本文件，画出下列分段函数所表示的曲面。（1）编写M脚本文件的步骤图 8.1-1-1 MATLAB Editor/Debugger 窗口,M文件
26、,M文件有两种形式：一是命令文件；另一种是M函数文件。,点击MATLAB指令窗工具条上的New File图标，就可打开如图8.1-1-1所示的MATLAB文件编辑调试器MATLAB Editor/Debugger。其窗口名为untitled，用户即可在空白窗口中编写程序。比如输入如下一段程序zx81.m%zx81.mThis is my first example.a=2;b=2;%clf;x=-a:0.2:a;y=-b:0.2:b;for i=1:length(y)for j=1:length(x)if x(j)+y(i)1 z(i,j)=0.5457*exp(-0.75*y(i)2-3.7
27、5*x(j)2-1.5*x(j);elseif x(j)+y(i)=-1 z(i,j)=0.5457*exp(-0.75*y(i)2-3.75*x(j)2+1.5*x(j);,M脚本文件,else z(i,j)=0.7575*exp(-y(i)2-6.*x(j)2);end endendaxis(-a,a,-b,b,min(min(z),max(max(z);colormap(flipud(winter);surf(x,y,z);点击编辑调试器工具条图标，在弹出的Windows标准风格的“保存为”对话框中，选择保存文件夹，键入新编文件名（如zx81），点动【保存】键，就完成了文件保存。（2）运
28、行文件使zx81.m所在目录成为当前目录，或让该目录处在MATLAB的搜索路径上然后运行以下指令，便可得到图形。zx81,M脚本文件,Matlab zx81,M脚本文件,M函数文件,M函数是由function语句引导的，其基本结构如下：,1）function返回变量列表函数名(输入变量列表)注视说明语句段，由引导输入、返回变量格式的检测函数体语句,2）还可以运用nline()函数来直接编写程序，其调用格式为：fun=inline(函数内容，自变量列表),例如：,可以用f=inline(sin(x.2+y.2),x,y)直接定义。,3）还可以借助匿名函数，其定义形式更为简洁,调用格式如下：f=(
29、变量列表)函数内容,例如：f=(x,y)sin(x.2+y.2),练习：应用上面给出的三个定义函数的方法计算函数,在(2,3)处的函数值,M函数文件,设可逆方阵A，试编写同时求 的M函数文件。,在编辑窗口输入以下的代码,Function da,a2,inva,traa=comp4(x)da=det(x);a2=x2;inva=inv(x);traa=x,保存该文件，文件名为comp4.M,在命令窗口输入下列语句：,M文件的调用,A1 2;5 8;comp4(A),Matlab comp4,M函数文件,题目1：Fibonacci数组的元素满足Fibonacci 规则：，；且。现要求该数组中第一个
30、大于10000的元素。用for循环指令来寻求Fibonacc数组中第一个大于10000的元素。题目2：利用while循环来计算1!+2!+50!,两个编程的小例子,a(1)=1;a(2)=1;i=2;while a(i)=10000 a(i+1)=a(i-1)+a(i);%当现有的元素仍小于10000时，求解下一个元素。i=i+1;end;i,a(i),i=21 ans=10946,两个编程的小例子,n=100;a=ones(1,n);for i=3:n a(i)=a(i-1)+a(i-2);if a(i)=10000 a(i),break;%跳出所在的一级循环。end;end,i ans=1
31、0946i=21,两个编程的小例子,题2编程：sum=0;i=1;while i51prd=1;j=1;while j=i prd=prd*j;j=j+1;endsum=sum+prd;enddisp(1!+2!+50!的和为:)sum,两个编程的小例子,多项式拟合,例如：对下表数据进行多项式拟合,x=1 2 3 4 5 6 7 8 9y=9 7 6 3-1 2 5 7 20;p=polyfit(x,y,3);xi=0:0.2:10;yi=polyval(p,xi);plot(xi,yi,x,y,r*),函数插值,1）有一个 矩阵,编程求出其最大值及其所处的位置.,今天的作业,2）编程求,3）有一函数,写一程序,输入自变量的值,输出函数值.,4)在word文档中直接计算完成：)生成向量a=1 2 3 4 5 6 7 8 9 8 7 6 5 4 3 2 12)生成矩阵b,b的每行元素就是向量a，矩阵b为17行17列的方阵3）将b的第三行元素去掉,
