数学实验与数学建模（matlab在建模中的应用）

第六章 数学实验与数学建模 学习目标 1.掌握利用Matlab软件进行了相关的数学运算的方法. 2.以软件辅助来完成数学实验. 3.了解数学建模思想方法，能够对一些简单问题建立数学模型求解分析. 教学要求 能力模块 能力要求 相关 知识点 高中化学知识点免费下载体育概论知识点下载名人传知识点免费下载线性代数知识点汇总下载高中化学知识点免费下载 运算能力 要求学生知道数学中运算所对应的相关 (1)数学的相关知识 Matlab基本函数及其用法，并能够运用相 (2)Matlab软件的语法 关函数完成数学的基本运算。 (3)相关函数的用法 实验能力 要求学生能够以软件作为辅助工具，按 线性代数中的相关数学理 要求完成相关实验 论与思想方法 建模能力 要求学生了解数学建模的思想方法，具备 一些基本的数学建模的 一定的数学建模能力 方法 Matlab是Mathworks公司推出的用于数值计算的交互式软件系统，具有强大的数值分析、矩阵运算、信号处理、图形显示和建模仿真功能. Matlab是“Matrix Laboratory”的缩写，意思是“矩阵 实验室 17025实验室iso17025实验室认可实验室检查项目微生物实验室标识重点实验室计划 ”，其强大的数据处理能力和丰富的工具箱使它的编程极为简单， 因此，它成为科学家和 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 技术人员解决实际问题的首选计算工具软件。 本章的第一节主要介绍Matlab软件的简单使用方法，从第二节到第六节在讲解Matlab用于解决高等数学和线性代数中的相关计算的函数基础上, 通过一些简单的数学实验例题,让学生体会如何用Matlab辅助解决数学问题. 最后，通过一些与线性代数相关的数学建模实例，让学生掌握数学建模的简单方法,学会利用Matlab软件辅助解决实际问题，以培养学生良好的数学意识和数学素质. 6.1 Matlab环境及使用方法 6.1.1 Matlab窗口管理 Matlab启动后显示三个窗口，如图6.1所示。左上窗口为工作区间窗口，显示用户定义的变量及其属性类型及变量长度。工作区间窗口也可显示为当前目录窗口，显示Matlab所使用的当前目录及该目录下的全部文件名。左下窗口为历史窗口，显示每个工作周期(指Matlab启动至退出的工作时间间隔)在命令窗口输入的全部命令，这些命令还可重新获取应用。右侧窗口为Matlab命令窗口，可在里面输入相关运算命令，完成相应计算。三个窗口中的记录除非通过菜单下的清除操作，否则将一直保存。 Edit 1 图6.1 Matlab运行期间(即程序退出之前),除非调用Clear函数，否则Matlab会在内存中保存全部变量值，包括命令输入的变量以及执行程序文件所引入的变量。清除工作空间变量值也可以通过下拉菜单中的Clear Workspace命令实现。Clear函数可以清除内存中Edit 的所有变量。 Matlab命令窗口输入的信息会保持在窗口中，并可通过通过滚动条重新访问。一旦信息量超出其滚动内容容量，则最早输入的信息将会丢失。可以通过在命令窗口中输入clc命令来清除命令窗口中的内容，也可以通过下拉菜单中的Clear Command Window子菜Edit 单清除，但这个操作仅清除命令窗口中的内容，但不能删除变量，要删除变量，只能通过clear. 为在命令窗口中能够更加清晰地显示字母及数字，Matlab提供了format函数的几种功能。其中 format short e 是系统默认的设置，显示5位数字， format long e 显示格式转换为16位数字加3位指数位的显示格式. format compact preferences命令将剔除显示中多余的空行或空格. 这此属性值也可通过单击菜单的File子菜单，弹出的preferences设置窗口后选择Command Window项进行设置(如图6.2). 图6.2 2 使用Matlab过程中有两个有用的组合键“^c”(Crtl+c)用于终止程序或函数的执行，也可用于退出暂停的程序或函数;“^p”(Ctrl+p)用于将最近键入的信息显示在Matlab命令窗口中，按Enter键可再次执行该命令，连续按两次“^p”,可调用上两次的输入信息，以此类推. 6.1.2 Matlab的基本语法 Matlab允许用户创建的变量名不可超过63个字符，多余部分将被忽略掉.变量名要求以大写或小写字母开头，后面跟大小写字母，数字或下划线。字符间不允许有空格. 变量名区分大小写，例变量名A1与a1表示不同变量. 此外，不能使用希腊字母，或者上下标字符作 ,为变量名，但可以拼写希腊字母，或在下标字符前加入下划线表示变量. 例如，可写为1lamda_1. Matlab在命令窗口运行时，要求首先在“>>”提示符后定义一个或多个变量，并进行赋值，然后表达式才能够使用变量。赋值运算符为“=”，输入变量名和等号后，按Enter表示结束.例如要实现a,2,则要在命令窗口中进行以下信息交互: >> a=2 用户输入 a = 系统响应 2 注:表达式后加分号(;)可省略系统响应信息的显示. Matlab允许在一行中输入多个表达式，表达式间以逗号或分号进行分隔，行尾以Enter键结束.用逗号分隔时系统会回显输入的值，如果用分号分隔表达式，不会输出响应信息. 例如按如下格式输入信息: >> a=2;b=2.5,c=3; 系统显示为: b = 2.5000 此时变量a和b的值不显示，但内存中存在. 标量的加、减、乘、除和幂运算分别用+,- ,*,/,^表示，默认的运算次序为:幂运算为最高，其次为乘除，最后为加减。同时在表达式中可用圆括号来确定运算次序。 c3,, 例6-1计算当a,2,b,3,c,6,时的值. t,,,1，2ab,, 输入: >> a=2;b=3;c=6;t=(3/(1+2*a*b))^c 结果: t = 1.5103e-004 Matlab软件中包含了大量的函数，表6.1、表6.2、表6.3、表6.4给出了一些常用的函数及其功能 说明 关于失联党员情况说明岗位说明总经理岗位说明书会计岗位说明书行政主管岗位说明书 . 3 表6.1 常用初等函数的Matlab表示 数学函数 Matlab表示 exp(x) x e sqrt(x) x log(x) lnx log10(x) lgx abs(x) x sign(x) sgn(x) 表6.2 三角函数和双曲函数的Matlab表示 三角函数 双曲函数 数学函数 原函数 反函数 原函数 反函数 正弦 sin(x) asin(x) sinh(x) asinh(x) 余弦 cos(x) acos(x) cosh(x) acosh(x) 正切 tan(x) atan(x) tanh(x) atanh(x) 余切 cot(x) acot(x) coth(x) acoth(x) 正割 sec(x) asec(x) sech(x) asech(x) 余割 csc(x) acsc(x) csch(x) acsch(x) 表6.3 Matlab复数运算函数 Matlab函数 功能 Complex(a,b) 表示复数 a，bi zabs(z) 取复数的模 zconj(z) 取复数的共扼复数 zreal(z) 取复数的实部 zimag(z) 取复数的虚部 表6.4 小数转换为整数的Matlab函数 Matlab函数 功能 fix(x) 舍去小数部分 round(x) 四舍五入 ceil(x) 取上整 floor(x) 取下整 4 Matlab软件中还有一些内置的常量，如表6.5所示 表6.5 部分常量的Matlab表示 数学常量 Matlab表示 , pi , inf 虚数单位 i(或j) i 一些特殊字符在Matlab中有特殊的功能，如表6.6所示 表6.6 特殊字符及其功能说明 符号 名称 功能 (a)小数点 . 句号 (b)向量或矩阵的一种操作类型。例如 c,a.*b (a)参数分隔符 ， 逗号 (b)几个表达式在同一行时放在每个表达式之后 (a)放在表达式末尾不显示计算结果 ; 分号 (b)在创建矩阵的语句中指示一行的结束，例如: ,,m,xyz;abc : 冒号 (a)创建向量表达式分隔符，例如: x,a:b:c (b)对矩阵而言，A(:,k)表示第列所有元素; kA A(k,:)表示第行所有元素 k ( ) 圆括号 z(a)矩阵中某一元素的下标指示，如z(j,k)表示矩阵行列的元素 jk (b)算术表达式分隔符，如a^(b，c) (c)函数参数分隔符，如sin(x) 创建一组数值、向量、矩阵或字符串(字母型) [ ] 方括号 创建单元矩阵或结构 { } 大括号 注释分隔符. 用于指示注释的开始，MATLAB编译器会忽略其右边的内容。 % 百分号 但用于一对引号内部定义字符串时除外，如:a='pl=14% of the totle' (a)'Expression'表明Expression为字符串(字母型) ‘ 引号 (b)表示向量或矩阵的转置 作为数据创建语句的分隔符，如; ,,c,ab 空格 或者作为字符串语句的一个字符 有了上述基本知识后，我们就可以利用Matlab进行一些简单的运算了. ,2sin，arcsin1,eln7例6-2 计算 3 输入: >> sin(pi/3)+asin(1)-exp(2)*log(7) 结果: ans = -11.9416 5 注:命令窗口作为计算器应用且未将计算结果分配给表达式时，Matlab默认将计算结果分配给变量名ans. 6.1.3 Matlab的编辑器 为了避免在命令窗口中输入多个表达式或重复输入一组表达式，且输入的表达式在执行完后能够进行修改， Matlab定义了程序文件:包含一组命令的文件，文件中每一命令行的执行都同在命令窗口里的执行一样.程序文件可由Matlab编辑器创建，编辑器的启动可由 下拉菜单中子菜单下的来完成，窗口如图6.3所示. M,fileFileNew 图6.3 使用Matlab编辑器创建的程序文件，点击图6.3所示的Run图标，可执行文件. 但在执行前首先要保存，保存可用Saveas下拉菜单中的子菜单，文件扩展名“.m”将自File 动添加到文件名末尾. 文件名的命名规则与变量名的命名规则相同. 例6-3 编写一个m文件，用来计算下面的表达式 L,2Scos,，,(R，r)，2,(R,r) R,r其中R,30,r,12,S,50,,arcsin()，，求 LS m文件如下: R=30; r=12; S=50; theta=asin((R-r)/S); L=2*S*cos(theta)+pi*(R+r)+2*theta*(R-r) 运算结果: >> L = 238.4998 Matlab编辑器的另外一个重要的功能是用来定义函数. 函数具有模块化、结构化和重用性的特点，是Matlab程序的重要组成部分，下面我们给出Matlab定义函数的方法. Matlab定义函数格式如下: Function [OutputVariables]=FunctionName(InputVariables) 6 %Comments Expression(s) 其中，OutputVariables是输出变量，变量之间用逗号隔开;InputVariables是输入变量，变量之间用逗号隔开;Function是保留字;FunctionName是函数名称，命名规则与变量命名规则一致;函数文件可以存储在任何已定义或即将定义的路径下，但文件名必须是FunctionName.m; %Comments是函数程序的注释部分，为了增强程序的可读性;Expression(s)是实现函数功能的主体程序，由一行或多行Matlab语句组合完成，函数的定义是在Matlab的编辑器中完成的. 注:当函数只是为了实现某一特定功能，而不涉及参数传递时，OutputVariables或InputVariables是可以省略的. ,2f(x,y),x，sinxy例6-4 定义函数，并求当时的函数值. x,5,y,2 定义函数程序如下(在编辑器中输入): function value1=hanshu1(x,y) value1=x^2+sin(x*y); 保存为hanshu1.m 在命令窗口中输入如下信息: >> hanshu1(5,pi/2) 结果: ans = 26 6.2 高等数学与Matlab Matlab的符号运算工具箱包含了微积分运算、解微分方程等几个方面的工具，本节我们主要介绍利用Matlab软件辅助解决高等数学中的相关计算问题. 6.2.1 求函数极限 求函数极限是高等数学中的重要计算，它在数学模型的分析和预测中起着重要的作用，Matlab提供了求极限的符号函数limit来完成这项工作，其调用格式和功能如表6.7所示. 表6.7 limit函数使用说明 格式 功能 x,al=limit(S,x,a) 求符号函数，当时的极限 S l=limit(S,x,a,’right’) ，x,a求符号函数，当时的右极限 S l=limit(S,x,a,’left’) ,x,a求符号函数，当时的左极限 S 注:对于符号函数S，也可以是普通的Matlab内置或自定义函数. 下面通过实例来说明这一函数用法: xxlime和lime例6-5 求 ,，,,,,xx 输入: >> syms x 7 >> r_value=limit(exp(x),x,+inf), l_value=limit(exp(x),x,-inf) 结果: r_value = inf l_value = 0 注:“syms x”语句是用来创建符号变量. 1xxlime和lime例6-6 求 ,x,0x,0 输入: >> syms x >> l_value=limit(exp(1/x),x,0,'left'),value=limit(exp(1/x),x,0) 结果: l_value = 0 value = NaN 注:NaN在Matlab里面表示数值(结果)不存在. 6.2.2 求导运算 求函数的导数和微分是高等数学中的基本计算，也是解决实际问题过程中的常用的一种 方法，Matlab提供了专门求导数的符号函数diff，能够帮助解决求导数和微分的问题.diff 的调用格式和功能如表6.8所示 表6.8 diff函数使用说明 格式 功能 xd=diff(S,’x’) 求符号函数关于变量的一阶导数(偏导数) S xnd=diff(S,’x’,n) 求符号函数关于变量的阶导数(偏导数) S 注:对于符号函数S，也可以是普通的Matlab内置或自定义函数. 下面我们通过两个实例来说明这函数的用法: 3,sin(xy),sin(xy)d,例6-7 求导数sinx，求偏导数 3,x,xdx 输入: >> df1=diff('sin(x)','x'),df2=diff('sin(x*y)','x'),... df3=diff('sin(x*y)','x',3) 结果: df1 = cos(x) df2 = cos(x*y)*y df3 = -cos(x*y)*y^3 8 或输入: >> syms x y >> df1=diff(sin(x),x),df2=diff(sin(x*y),x),... df3=diff(sin(x*y),x,3) 也可得到相同的结果，后面类似问题两种输入方法均可，不再赘述. 2xf(x,y),x，sinxy例6-8 求函数对变量的一阶、二阶偏导. 在例6-4定义的函数基础上，我们在命令窗口中输入如下信息: >> syms x y >> df1=diff(hanshu1(x,y),x),df2=diff(hanshu1(x,y),x,2) 结果: df1 = 2*x+cos(x*y)*y df2 = 2-sin(x*y)*y^2 多元函数的雅克比矩阵是一个非常重要的概念，多元函数的积分、矩阵微积分及坐标转 换等都要用到它，在Matlab中用jacobian函数求雅克比矩阵，其调用格式如下: j=jacobian(F,V) 其中是函数向量，是自变量向量. VF 2tTf,[xsin(xyz)e]例6-9 计算函数向量关于自变量的雅克比矩阵. x,y,z,t 在命令窗口中输入下列命令: >> syms x y z t; >> f=[x^2 sin(x*y*z) exp(t)]; >> v=[x y z t]; >> j=jacobian(f,v) 输出结果为: j = [ 2*x, 0, 0, 0] [ cos(x*y*z)*y*z, cos(x*y*z)*x*z, cos(x*y*z)*x*y, 0] [ 0, 0, 0, exp(t)] 注:关于向量、矩阵的表示及相关运算，我们将在下一节作详细介绍，这里主要介绍一 下jacobian函数的用法. 6.2.3 积分运算 Matlab提供了两类计算积分的函数，一类是符号积分函数，一类是数值积分函数.下面 我们先介绍符号积分函数int，其调用格式和功能如表6.9所示: 表6.9 int函数使用说明 格式 功能 int(F) 求符号表达式求关于自变量的不定积分 F sint(F,’s’) 求符号表达式求关于自变量的不定积分 F aint(F,a,b) 求符号表达式当自变量从到的定积分 bF saint(F,’s’,a,b) 求符号表达式当自变量从到的定积分 bF 9 下面我们通过实例来说明函数的用法: 211,x例6-10 计算和 f1,ysinxdxf2,dxxydy,,,00 输入: >> syms x y >> f1=int(y*sin(x),x),f2= int(int(x*y,y,0,sqrt(1-x^2)),x,0,1) 结果: f1 = -y*cos(x) f2 = 1/8 下面我们来介绍数值积分函数quad,dblquad,triplequad的用法 表6.10 数值积分函数使用说明 格式 功能 aquad(F,a,b) 从到计算符号函数的数值积分 bF dblquad(F,x1,x2,y1,y2) ,,x,x,y,y计算符号函数在区域F1212 上的二重数值积分 triplequad(F,x1,x2,y1,y2,z1,z2) 计算符号函数在区域F ,,x,x,y,y,z,z上的三重数值积分 121212 注:关于上述三个函数更详细的用法请参考软件提供的帮助. 1112,f1,dxysinxdyf2,dxdy(ysinx，zcosx)dz例6-11 计算和 ,,,,,,00,110 输入: >>f1=dblquad('sqrt(y)*sin(x)',-1,1,0,2),... f2=triplequad('y*sin(x)+z*cos(x)',0,pi,0,1,-1,1) 结果: f1 = -4.4306e-018 f2 = 2.0000 6.2.4 微分方程的求解 Matlab提供了两类求解微分方程的函数，一类是求微分方程的符号解，一类是求微分 方程的数值.本部分主要介绍利用dsolve求微分方程的符号解，其调用格式如表6.11. 表6.11 dsolve函数使用说明 格式 功能 dsolve(’equation’) 求方程equation通解，自变量默认为 t dsolve(’equation’,’condition’) 求方程equation满足初始条件condition 的特解，自变量默认为 t xdsolve(’equation’,’x’) 求方程equation通解，自变量为 10 dsolve(’equation’,’condition’,’x’) 求方程equation满足初始条件condition x的特解，自变量为 n"f(x,y,y",？y),0y说明:(1)equation代表常微分方程式，且用Dy表示，用y'D2y表示，以此类推. Condition为初始条件. (2)函数dsolve可用来求常微分方程、方程组的解. 3yx，',y例6-12 求方程的通解. x 输入: >> dsolve('Dy=(y+x^3)/x') 结果: ans = -x^3+exp(1/x*t)*C1 x这显然是个错误的结果，程序将看作是自变量，而是常量了. t 正确输入如下: >>dsolve('Dy=(y+x^3)/x','x') 结果: ans = 1/2*x^3+x*C1 2xy",2y',3y,3xe例6-13 求方程满足的特解. y,0,y',1x,0x,0 输入: >>solve('D2y-2*Dy-3*y=3*x*exp(2*x)','y(0)=0,Dy(0)=1','x') 结果: ans = -1/3*exp(2*x)*(3*x+2)-1/3*exp(-x)+exp(3*x) 求常微分方程数值解方面，Matlab提供了比较丰富的函数,其中常用的为ode45和ode32,两个函数都运用了基本的龙格-库塔(Runge-Kutta)数值积分法的变形,其基本调用格式如下: [t,Y]=ode45(diffeq,[t0 tn],y0,options,arg1,arg2,...) 其中，第一个参数diffeq是一个m文件的名字，该文件用来计算微分议程的右端表达式;t0,tn定义了求解总区间，如果t0未知，则假设其为0;y0为初始条件;options为函数ode45参数;而参数arg1,arg2，...是传给用户自己定义的diffeq程序的. ode32用法与ode45相同. 求常微分方程的数值解问题现阶段我们遇到不是很多，所以在这里不作详细讲解，如果用到可查阅Matlab帮助文件. 6.3 线性代数与Matlab 前面我们的例子中的变量都是标量，而Matlab的强大之处在于它是以数组作为变量来进行运算的，这将使我们处理大量的数据变得十分方便，本节我们将重点介绍如何定义和使用Matlab中的矩阵变量来进行简单的计算. 6.3.1 向量的创建及对元素的操作 a,[a1,a2,？,an]或[a1a2？an]行向量在Matlab表示为:,列向量在Matlab中 11 表示为b,[b1;b2;？;bn]，其中的可以是数值、字符串、预先定义的变量或表达式. ai或bi 向量的创建可以是直接按格式输入元素，也可以用Matlab函数完成，这里面有两个常 用的创建向量的命令分别是冒号(:)和linspace，其调用格式如下 (1) x,x0:d:x1 n其中表示向量的起始值，表示步长，表示向可能结束值，若向量的元素为个，x0dx1 则它们满足 . x0，nd,x1 x (2) =linspace() x0,x1,n nn 它表示在区间[x0,x1]上等距离的取个数，作为向量的值， 系统默认为100. 例如输入如下信息: >>x1=0.2:2:12,x2=linspace(0.2,12,6) 结果: x1 = 0.2000 2.2000 4.2000 6.2000 8.2000 10.2000 x2 = 0.2000 2.5600 4.9200 7.2800 9.6400 12.0000 对于向量及其元素的操作见表6.12所示 表6.12 向量及其元素操作 格式 功能 n=length(x) 返回向量元素个数 xx’ 求向量的转置 xx(i) 取向量的第个元素 i xax(i)=a 取向量的第个元素赋值为 i xx(end) 取向量的最后一个元素 y=x([i,j,k]) x取的第,元素构成新的向量 j,ki x[xnew,index]=sort(x) xnew是升序排列后的向量，index是 xxnew中元素在中的位置构成的向量 xindex=find(condition) 将向量中满足condition条件的元素 下标做成向量index x[xmax,kmax]=max(x) xmax表示向量的最大值，kmax表示 最大值位置 [xmin,kmin]=min(x) 功能与max类似 f(x) 对向量中每个元素求函数值 xsum(x) 求向量的所有元素之和 xcumsum(x) 返回以向量的元素的部分和数列为元 素的向量 norm(x) 求向量的范数 mean(x) 求向量的平均值 std(x) 求向量的 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 差 在对向量及其元素的操作过程中，上述命令总是组合使用. 12 例6-14 在区间[0,2,]上等间距取8个值，做成行向量，求 (1)求各点处正弦值，取出其中最大值，并指明位置 (2)取出所有大于0的正弦值，将其加上1后，做成一个新的列向量. 为一更好的完成上述工作，我们将其做成一个m文件，程序如下: x=linspace(0,2*pi,8) %取向量x y=sin(x) %对x中的每一个元素计算sin值,新向量赋值给向量y [ymax,kmax]=max(y) %求y中的最大值 index=find(y>0) %将y中大于0的元素下标做成一个向量 y1=y(index)+1 %利用下标向量index将y中小于0的元素取出加1 y1=y1’ %将行向量转置为列向量，重新赋值给y1 保存程序后运行，结果为: x = 0 0.8976 1.7952 2.6928 3.5904 4.4880 5.3856 6.2832 y = 0 0.7818 0.9749 0.4339 -0.4339 -0.9749 -0.7818 -0.0000 ymax = 0.9749 kmax = 3 index = 2 3 4 y1 = 1.7818 1.9749 1.4339 y1 = 1.7818 1.9749 1.4339 6.3.2 矩阵的创建及对元素的操作 一个3行4列的矩阵在Matlab表示为 A =[a11,a12,a13,a14;a21,a22,a23,a24; a31,a32,a33,a34;a41,a42,a43,a44] A 或 =[a11,a12,a13,a14;„ A a21,a22,a23,a24;„ a31,a32,a33,a34;„ a41,a42,a43,a44] a其中可以是数值、字符串、预先定义的变量或表达式,上述格式中的逗号可以用空格代替.ij 简而言之，矩阵在Matlab中表示为用[ ]括起来的一组元素，同一行元素之间用逗号(,)或空隔分隔，行与行之间用分号(;)分隔. 矩阵的创建可以是直接按格式输入元素，也可以用Matlab函数完成，表6.13 给出了创建特殊矩阵的函数 13 表6.13 创建特殊矩阵函数 格式 功能 m，nones(m,n) 创建阶所有元素都是1的矩阵 m，nzeros(m,n) 创建阶所有元素都是0的矩阵 nadiag(a) 创建阶对角阵，对角线上元素为向量的 元素 neye(n) 创建阶单位矩阵 nmagic(n) 创建阶魔方阵 m，nxrepmat(x,m,n) 创建阶所有元素都是的矩阵 [U,V]=meshgrid(s,t) 由两个行向量，构建两个维数相同的s和t 矩阵和 UV x注:repmat参数中的可以是标量、向量或矩阵. 为了能更好的说明meshgrid的功能，我们给出下面一个简单的例子: 如果，则命令:[X,Y]=meshgrid(x,y)将返回的两x,[x1,x2,x3,x4],y,[y1,y2,y3] 个的矩阵: 3，4 y1y1y1y1x1x2x3x4，，，， ,,,,Y,y2y2y2y2X,x1x2x3x4 ， ,,,, ,,,,x1x2x3x4y3y3y3y3，，，， 这个函数在做三维图形时会有十分有用. 向量可以看作是列(行)矩阵，而矩阵的每一列(行)也可以看作是一个向量，所以我 们也可以通过向量的组合来创建矩阵,比如下面的命令就可以创建一个新的矩阵: 输入: >> A=[2:3:15;linspace(10,21,5);ones(1,5)] 结果: A = 2.0000 5.0000 8.0000 11.0000 14.0000 10.0000 12.7500 15.5000 18.2500 21.0000 1.0000 1.0000 1.0000 1.0000 1.0000 对于矩阵及其元素的操作见表6.14所示 表6.14 矩阵及其元素操作 格式 功能 a取出矩阵的第行与第列交叉点元素 jiAA(i,j) ij a,a赋值操作 A(i,j)=a ij 取出矩阵的第行所有元素，构成一个行向iAA(i,: ) 量 取出矩阵的第列所有元素，构成一个列向jAA( :,j) 量 14 的行和行，与列和列和取出矩阵j1j2i1i2A A([i1,i2],[j1,j2,j3]) 列交叉点的元素，构成的一个子矩阵 j3A sort(A) 按升序排列矩阵列元素 max(A) 求矩阵中每一列的最大值 A min(A) 求矩阵中每一列的最小值 A 将矩阵的各列元素相加，返回一个长度等于Asum(A) 列数的行向量 A mn[m,n]=size(A) 返回矩阵的行数和列数 A f(A) 对的每个元素求函数值 A 注:A(i,,=[]或A(:,j)=[]，则表示删除第行或第列. ji 1234，， ,,A,5678例6-15 给定矩阵,完成如下工作 ,, ,,9101112，， (1)将矩阵的第4列全部赋值为0元素. A (2)求矩阵的最大值. A (3)取A的第1，3行和2，3列交叉点元素，取自然对数后，构成新的矩阵. B我们将其做成一个m文件，程序如下: A=[1,2,3,4;5,6,7,6;9,10,11,12]; %创建矩阵A A(:,4)=0 %A的第4列全部赋值为0 maxvalue=max(max(A)) %求矩阵A的最大值 B1=A([1,3],[2,3]) %B1为A的第1，3行，第2，4列元素构成的矩阵 B=log(B1) %B中的每个元素取自然对数 保存后，运行结果如下: A = 1 2 3 0 5 6 7 0 9 10 11 0 maxvalue = 11 B1 = 2 3 10 11 B = 0.6931 1.0986 2.3026 2.3979 6.3.3 矩阵的运算 矩阵的运算主要分成两类，矩阵的点运算和线性代数中的矩阵运算. 同时需要指明的 是，因为向量可以看成是特殊的矩阵，所以本部分所讲的关于矩阵的各种运算，同样适用于 向量. 下面我们分别将两类运算做一介绍. 在Matlab中，矩阵的点(.)运算是指对同阶矩阵中逐个元素进行的算术运算，分别为点 15 乘(.*)、点除(./)和点幂(.^).下面我们以两个矩阵来分别说明其功能. 2，3 xxxyyy，，，，111213111213X令，同时Y ,,,,,,xxxyyy212223212223，，，， 则 x*yx*yx*y，，111112121313Z1X.*Y ,,,,x*yx*yx*y212122222323，， x/yx/yx/y，，111112121313Z2X./Y ,,,,x/yx/yx/y212122222323，， x^yx^yx^y，，111112121313Z3X.^Y ,,,,x^yx^yx^y212122222323，， X,xxZ1,x*Y这里需要说明的是对于点乘，如果且为标量，则点可以省略，即，000 Y,yyZ1,X*yY,yy类似的，如果且为标量，则;对于点除，如果且为标量，00000 Z2,X/yX,xx则点可以省略，即，但是如果且为标量，则点不可以省略，000Z2,x./YX,xY,y;对于点幂运算，不论是还是取标量，点都不可以省略，即000Z3,x.^YZ3,X.^y和 . 00 x例6-16 求在区间区间内，等距离取6个值时，表达式 0,x,2 bx，ccos(),axy,e x，c ,0.3,1.2,的值. 其中，. a,b,c,3 M文件的程序代码如下: a=0.3;b=1.2;c=pi/3; x=linspace(0,2,6); y=exp(-a*x).*cos(b*x+c)./(x+c) 结果: y = 0.4775 0.0267 -0.1800 -0.2463 -0.2302 -0.1718 例6-17 利用向量元素操作函数中的sum和cumsum，求表达式 5nz,n ,n,1 的值，及其部分和数列对应向量. 输入: >> n=1:5;„ an=n.^n,„%计算级数一般项 16 s=sum(an),„%求表达式值 sn=cumsum(an)%求部分和 结果: an = 1 4 27 256 3125 s = 3413 sn = 1 5 32 288 3413 在线性代数中矩阵的运算包括:加、减、乘、转置、求逆等，其基础理论在前面都做详 细的介绍，这里不再赘述。Matlab对于矩阵的这些运算都能够实现，表6.15给出相关函数 及其使用说明. 表6.15 矩阵数学运算的函数 格式 功能 C=AB C表示矩阵与的和(差) ,AB C=A*B C表示矩阵左乘矩阵 AB 求方阵的行列式 AAdet(A) inv(A)或A^(-1) 求矩阵的逆矩阵 A lamda返回的是矩阵的特征值，是特征值AP[P,lamda]=eig(A) 对应的单位特征向量构成的正交矩阵 rank(A) 确定矩阵线性无关的行或列的值(矩阵的秩) A’ 返回矩阵的转置 A A^k 返回矩阵的次幂 kA trace(A) 返回矩阵的迹 A poly(A) 返回的是矩阵的特征多项式(系数向量) A rref(A) 将化成阶梯形行最简形式(与等价) AA orth(A) 将非奇异矩阵正交化规范化 A flipud(A) 将矩阵按列旋转180º A fliplr(A) 将矩阵按行旋转180º A 注:如果矩阵和是行向量和列向量，则表示内积，也可以用C=dot(A,B)C,A*BAB 表示. 232f(x,x,x),x，2x，5x，2xx，6xx，2xx例6-18 化二次型用正交变换化123123121323 为标准形，并写出变换矩阵. 113，， ,,A,121因为二次型矩阵，所以先用Matlab求特征值和正交矩阵如下: ,, ,,315，， >> A=[1 1 3;1 2 1;3 1 5];... [Q,lamda]=eig(A) 结果: Q = 0.8835 -0.0253 0.4677 17 -0.1667 -0.9501 0.2636 -0.4377 0.3108 0.8437 lamda = -0.6749 0 0 0 1.6994 0 0 0 6.9754 222所以二次型的标准形为，变换矩阵为. Qf,,0.6749y，1.6994y，6.9754y123 求的逆矩阵为: Q >> inv(Q) 结果为: ans = 0.8835 -0.1667 -0.4377 -0.0253 -0.9501 0.3108 0.4677 0.2636 0.8437 ,1TQ,Q 进一步我们可以验证正交变换的一个性质:. 同时，我们可得到 y0.8835,0.1667,0.4377x，，，，，，11,,,,,,,1Ty,y,Qx,Qx,,0.0253,0.95010.3108x. 22,,,,,, ,,,,,,y0.46770.26360.8437x33，，，，，， 通过下面的命令，我们还可以验证正交矩阵将矩阵对角化为矩阵lamda: QA >> Q'*A*Q 其结果是矩阵lamda. 注:因为Matlab进行的是数值计算，所以给出的解为数值解，而我们在教材上计算得出的解基本上都是解析解，故两者之间有一些区别. 6.3.4 方程组求解 Matlab对于非齐次线性方程组的求解有如下两种方法: AX,b 法1:X=A\b 或 X=inv(A)*b 法2:C=[A,b] %写出增广矩阵 D=rref(C) %D的最后一列元素就是所求的解 xxx,,，,0,234,x，x，x，x,6,1234例6-19 求方程组的解. ,xxxx2，4，,2,,11234, ,3x，x,2x，2x,31234, 法1程序如下: 输入: >> A=[0 -1 -1 1;1 1 1 1;2 4 1 -2;3 1 -2 2];... b=[0;6;-1;3];... X=A\b 18 结果: X = 2.0000 -1.0000 3.0000 2.0000 法2程序如下: >> C=[A,b];D=rref(C) 结果: D = 1 0 0 0 2 0 1 0 0 -1 0 0 1 0 3 0 0 0 1 2 则的最后一列就是方程组的解，与法1结果相同. D 对于齐次线性方程组,Matlab中的函数null可以用来求方程组的基础解系.其AX,0 调用格式如下: beta=null(A,’r’) %beta的列向量是方程组的有理基. AX,0注:没有参数r，返回的是小数形式表示的基，有参数r则小数用最接近的分数表示 x，x,x,x,0,1234,2x,5x，3x，2x,0例6-20 求方程组的基础解系. ,1234 ,7x,7x，3x，x,01234, 输入: >> A=[1 1 -1 -1;2 -5 3 2;7 -7 3 1];... format rat,... %指定有理格式输出 beta=null(A,'r') 结果: beta = 2/7 3/7 5/7 4/7 1 0 0 1 所以解空间的有理基为向量beta的列向量，进可以写出通解: >> syms c1 c2 >> x=c1*beta(:,1)+c2*beta(:,2) 结果: x = [ 2/7*c1+3/7*c2] [ 5/7*c1+4/7*c2] [ c1] [ c2] 6.4 多项式的运算 19 多项式在Matlab中表示为用降幂排列的多项式的系数构成的向量，比如多项式42P,x，3x，x,6在Matlab表示为: >> P=[1,0,3,1,-6]; 如果将多项式的向量表示形式转变为符号表示形式，则采用函数ploy2sym, 比如将多 项式P在Matlab在中表示为符号形式，则 输入: >> poly2sym(P,’x’) %后面的参数x，表示以x为自变量 结果: ans = x^4+3*x^2+x-6 多项式的运算在Matlab中都表示为以系数向量为参数的运算，下面我们在讨论多项式 时都以向量指代，即系数向量P就表示多项式P. 多项式的常用运算见表6.16 表6.16 多项式运算 格式 功能 与和，向量与维数必求多项式SSPPPS ,需相同，缺项补0 求多项式和的乘积，也是向量和SPPconv(P,S) 的卷积 S r求多项式除以的商和余项，也SqP[q,r]=deconv(P,S) 是向量的解卷积运算 求方程p(x),0的解 roots(P) polyder(P) 求多项式的一阶导数 P polyder(P,S) 求多项式与的乘积的一阶导数 SP 求多项式与的商的一阶导数，SP/SP [q,d]=polyder(P,S) 其中表示分子，表示分母 dq x,ay=polyval(P,a) 求当时，多项式的的值. P 将多项式运算与前面的矩阵与向量运算结合，就可以得到利用Matlab求矩阵特征值和 特征向量的方法，其步骤如下: (1) 用P=poly(A)求方阵的特征多项式系数向量; AP (2) 用lamda=roots(P)求特征多项式的全部特征根lamda. P (3) 用函数beta=null(A-lamda*E)求出的基础解向量beta中的每一个列向量就是 lamda对应的一个特征向量. ,31,1，， ,,A,,75,1例6-21 求的特征值和特征向量. ,, ,,,66,2，， m文件程序如下: clc; clear; 20 A=[-3,1, -1;-7,5,-1;-6,6,-2]; %输入矩阵A P=poly(A) %求矩阵A的特征多项式(步骤1) lamda=roots(P) %求矩阵A的特征根(步骤2) b1=A-lamda(1)*eye(3); beta1=null(b1,'r') %求第一个特征值对应的特征向量 b2=A-lamda(2)*eye(3); %b2=rref(b2,1e-7) %为保证矩阵b2的奇异性而设的语句 beta2=null(b2,'r') %求第二个特征值对应的特征向量 b3=A-lamda(2)*eye(3); %b3=rref(b3,1e-7) %为保证矩阵b3的奇异性而设的语句 beta3=null(b3,'r') %求第三个特征值对应的特征向量 结果: P = 1.0000 -0.0000 -12.0000 -16.0000 lamda = 4.0000 -2.0000 -2.0000 beta1 = 0.0000 1.0000 1.0000 beta2 = Empty matrix: 3-by-0 beta3 = Empty matrix: 3-by-0 对于第二个特征值，之所以找不到特征向量是由于计算误差而引起的，输入:det(b2)， 得到的结果是:ans =1.7799e-014.计算误差使得本应为零的行列式成了一个很小的数，使 得b2X=0没有非零解.解决的方法是将b2化为阶梯形式，因为在rref函数中可以人为地设 ,71，10定容差，例如把容差设为，令b2=rref(b2,1e-7)，这样变换出的行阶梯形矩阵与原 矩阵b2是等价的，但它将会把一些小的数看成零，从而保证矩阵的奇异性. 所以，在上述 程序中的b2=A-lamda(2)*eye(3)与beta2=null(b2,'r')两条语句之间插入语句 Tb2=rref(b2,1e-7)，再执行程序，可得与特征值-2对应的特征向量为beta2=[1,1,0].对 于b3的处理原理相同. 6.5 绘制函数图像 Matlab中基本的二维绘图命令为: plot(xdata,ydata,symbol) y,f(x)其中xdata和ydata对应于的自变量和函数取值的向量， y,f(x)(xdata(i),ydata(i))表示点的坐标，此命令用于将各个点用实线连接，描绘出 21 的图像. symbol是用单引号(’’)引出的字符或符号，用于表示点(线)的颜色或形状特 征，系统默认是实线连接各点，symbol的部分可选项见表6.17 . 表6.17 symbol的部分可选项 线的类型 线的颜色 点的类型 符号 说明 符号 说明 符号 说明 - 实线 y 黄色 . 点 -- 虚线 r 红色 o 空心圆 : 点线 g 绿色 * 星号 -. 横点 b 蓝色 s 正方形 k 黑色 p 五角形 此外，plot(xdata1,ydata1,symbol1，xdata2,ydata2,symbol2，„)还可以在一个坐标 系中画条曲线. 例6-22 在一个坐标系中画出的图像，的用红实y,sinx,y,cosx,x,[0,2,]sinx cosx线表示，用绿横点线表示，并将点用蓝星号标出. (2,,1) m文件程序如下: x=linspace(0,2*pi); y1=sin(x); y2=cos(x); plot(x,y1,'r',x,y2,'g-.',x(end),cos(x(end)),'b*') 保存后，执行结果如图6.4 图6.4 为人提高图形的表达能力，Matlab还提供了给图形作注解的函数，见表6.18. 表6.18 给图形作注解的函数 注解函数 功能 x规定轴和轴的范围 yaxis grid 显示网格线 gtext 在对应鼠标输入处添加文本 legend 标明不同曲线的符号与线型的含义 (x,y)在点处添加文本 text 22 title 在图形顶部添加标题 xxlabel 用字符串标注轴 用字符串标注轴 yylabel 以上函数的具体用法，详见Matlab帮助文件. Matlab中绘制三维曲线的函数为: plot3(xdata,ydata, zdata ,symbol) 该命令是将二维图形函数plot的特性扩展到了三维空间，函数格式除了包括第三维信息外，基本用法与二维作图函数plot相同. Matlab中绘制三维曲面的图像有多种方法，其中两个最基本的是mesh和z,f(x,y) surf,它们带有相同的输入参数. mesh绘制的是网状图，surf绘制的是有阴影着色的网状图，其基本调用格式为: mesh(xdata,ydata,zdata)或surf(xdata,ydata,zdata). 此外，在进行数据分析时，常需要在一个画面上建立几个坐标第，每个坐标系显示不同 m，n的内容. Matlab中提供了窗口拆分函数subplot(m,n,p),表示将一个画面分成个图形区域，代表当前区域号，每个区域可以及对应一个独立的图. p 例6-23 在一个画面上作四个图，第一个显示曲线y,x^2,x,[,2,2];第二个显示曲 ,0.1t,x,ecost ,,0.1t22y,esintt,[0,6,]z,2,x,y,x,[,5,5],y,[,5,5]线;第三个显示的网格, ,z,t, 图，第四个显示上述曲面的有阴影着色的网格图. m文件程序如下: subplot(2,2,1) x=-2:0.05:2; y=x.^2; plot(x,y) %第一个图 subplot(2,2,2) t=linspace(0,6*pi); x=exp(-0.1*t).*cos(t); y=exp(-0.1*t).*sin(t); z=t; plot3(x,y,z) %第二个图 subplot(2,2,3) x=linspace(-5,5,25); y=linspace(-5,5,25); [X,Y]=meshgrid(x,y); %meshgrid的用法前面有介绍 Z=2-X.^2-Y.^2; mesh(X,Y,Z) %第三个图 subplot(2,2,4) surf(X,Y,Z) %第四个图 结果如图6.5所示 23 图6.5 关于Matlab提供的其它作图函数及修饰图形参数用法，如果需要，可查阅Matlab帮助文件，再这里不做过多讲解. 6.6 Matlab编程 Matlab程序实际上是将有关的Matlab命令为实现一定的功能而组合在一起，存储在一个文件中(m文件). 需要时，运行该文件，Matlab就会自动依次的执行文件中的命令，自动得到结果. 下面将编程所涉及的常用语句，包括分支语句、循环语句和转移控制语句做一简单介绍. 1. 分支---if语句 if语句为条件判断语句，通过满足不同的条件表达式将程序转向不同的部分. if语句调用格式如下: if 条件#1 表达式#1 elseif 条件#2 表达式#2 else 表达式3 end 其中“条件#n”通常应有一个(或多个)逻辑表达式，当表达式值为真时，执行“表达式#n”;否则跳过“表达式#n”. 在建立逻辑表达式时，常用的关系运算符如表6.19. 这里需要说明的是上述格式中的elseif和else语句是可选项，且程序中可同时出现多个elseif语句. 表6.19 常用关系运算符 符号 意义 == 等于 < 小于 > 大于 <= 小于等于 >= 大于等于 ~= 不等于 & 逻辑与 24 | 逻辑或 ~ 逻辑非 ,x,x,0, ,f(x),1,x,0例6-24 编写Matlab函数fun(x)，来表示分段函数，并计算, ,x,x,0, 的值. f(,2),f(0) 程序如下: function value=fun(x) %%定义分段函数 if x<0 value=-x; elseif x==0 value=1; else value=x; end 保存为fun.m后，在命令窗口中输入: >> y1=fun(-2),y2=fun(0) 结果: y1 = -2 y2 = 1 Matlab还提供了另外一个分支语句switch，用法比较简单,其能够实现的功能if语句 基本上都可实现，后面通过几个例子介绍了一下它的用法,这里不作详细讲解，如果需要深 入了解请查阅Matlab提供的帮助文件. 2. 循环---for语句 for循环可按指定次数重复执行一系列语句，其调用格式如下: for 循环变量=表达式1:表达式2:表达式3 循环体语句 end 其中“表达式1”的值为循环变量的初值，“表达式2”的值为步长，“表达式3”的值 为循环变量的终值，步长为1时，“表达式2”，可以省略. 例6-25 给定几位同学成绩chengji=[45,50,75,65,85,76,95,61,84,21,39,88 ],设计 一个分类器，输出为0表示不及格(<60)，1表示及格(>=60),-1表示分数错误. m文件程序如下: chengji=[45,50,75,65,85,76,95,61,84,21,39,88 ]; n=length(chengji); %求成绩向量长度 fenlei=zeros(1,n); %预先分配空间 for i=1:n if chengji(i)>=60&chengji(i)<=100 %判断 25 fenlei(i)=1; %分类 elseif chengji(i)<60&chengji(i)>=0 fenlei(i)=0; else fenlei(i)=-1; end end fenlei 执行结果: fenlei = 0 0 1 1 1 1 1 1 1 0 0 1 3. 循环---while语句 while循环无限次执行一条或多条语句，当条件不满足时停止循环。其调用格式如下: while 条件 循环体语句 end 其中，定义“条件”的表达式中可包含由“循环体语句”计算得出的一个或多个变量. SS的项数，表示如下: 例6-26 计算级数的部分和NNN N1S, ,N2nn,1 S求不超过的的最大值时，项数为多少, S,,/2NN m文件程序如下: format long %显示格式转为16位 S=pi/2; Sn=0;k=0; %给部分和Sn和循环变量k赋初值 while Sn> x=input('please enter x:') 结果: please enter x: 如果输入向量[1 2 3]，并按回车，则输出x的结果如下: x = 1 2 3 Matlab一般用disp命令和fprintf命令来显示结果. disp命令显示字符串时，字符串要用单引号(’’)括起来，并且用disp显示变量或字符串内容时，无“x=”或者“ans=”,例如继续在命令窗口中输入: >> disp(x) 结果: 1 2 3 fprintf命令能够以加后缀的方式显示不同的格式，fprint的格式说明符如表6.21所示: 表6.21 fprintf常用的类型说明符和特殊符 类型符 类型 字符 含义 %c 字符型 \n 换行 %s 字符串型 \t 水平制表(跳到下一个tab位置) %d 十进制整型 \b 退格，将当前位置移到前一列 %f 浮点型 \r 回车 %e 十进制指数型 \f 换页，将当前位置移到下页开头 %x 十六进制整数型 %% % %bx 十六进制浮点型 例6-27 在例6-25工作基础上，要求能够从键盘直接输入学生成绩向量，计算机直接输出学生及格与否. 27 m文件程序如下: clc,clear %清除命令窗口和内存变量 disp('这是一个学生成绩分类器,自动输出学生的分类情况') chengji=input('请输入学生成绩向量chengji='); n=length(chengji); %求成绩向量长度 fenlei=zeros(1,n); %预先分配空间 for i=1:n if chengji(i)>=60&chengji(i)<=100 %判断 fenlei(i)=1; %分类 elseif chengji(i)<60&chengji(i)>=0 fenlei(i)=0; else fenlei(i)=-1; end end disp('成绩结果为:') for i=1:n switch fenlei(i) %用switch分支结构判断 case 1 fprintf('第%d名学生成绩及格\n',i) %屏幕输出学生成绩判定结果 case 0 fprintf('第%d名学生成绩不及格\n',i) case -1 fprintf('第%d名学生成绩录入错误\n',i) end end 保存并运行程序，命令窗口中信息如下: 这是一个学生成绩分类器,自动输出学生的分类情况 请输入学生成绩向量chengji= 输入[45,50,75,65,85,76,95,61,84,21,39,188 ]，然后按回车键后，显示结果如下: 成绩结果为: 第1名学生成绩不及格 第2名学生成绩不及格 第3名学生成绩及格 第4名学生成绩及格 第5名学生成绩及格 第6名学生成绩及格 第7名学生成绩及格 第8名学生成绩及格 第9名学生成绩及格 第10名学生成绩不及格 第11名学生成绩不及格 第12名学生成绩录入错误 28 6.7 数学实验与数学建模 本节将从具体实例出发，通过详细讲解数学模型的建立和求解的过程，使学生掌握数学建模的思想方法，能够通过Matlab辅助计算，完成数学模型的求解分析过程. 为了与教材配套,本节数学模型所涉及的数学理论知识主要以线性代数为主,通过理论与实践相结合,来提高学生的数学应用能力，强化学生数学素质的培养. 29