matlab软件与数学建模.ppt

MATLAB软件与数学建模 主讲 鲜思东Email xiansdlearning 重庆邮电大学 常用数学建模软件 1 MATLAB软件 Mathematic软件2 LINDO LINGO软件3 SAS Spss Eviews Stata Splus R 统计分析 4 EXCEL软件5 其他 如CPLEX等 数学建模的基本方法和步骤 基本方法 根据对客观事物特性的认识 找出反映内部机理的数量规律 机理分析 机理分析没有统一的方法 主要通过实例研究 CaseStudies 来学习 以下建模主要指机理分析 测试分析 将研究对象看作 黑箱 通过对量测数据的统计分析 找出与数据拟合最好的模型 二者结合 机理分析建立模型结构 测试分析确定模型参数 数学建模的一般步骤 形成一个准比较清晰的 问题 了解实际背景 明确建模目的 搜集有关信息 掌握对象特征 数学建模的一般步骤 模型假设 在合理与简化之间作出折中 针对问题特点和建模目的 作出合理的 简化的假设 尽量采用简单新颖的数学工具 模型构成 用数学的语言 符号描述问题 发挥想象力 使用类比法 数学建模的一般步骤 模型求解 各种数学方法 数学软件和计算机技术 如结果的误差分析 模型对数据的稳定性分析 模型分析 与实际现象 数据比较 检验模型的合理性 适用性 模型检验 模型应用 Matlab与建模基本内容 一 MATLAB简介及基本运算二 MATLAB图形绘制功能三 MATLAB程序设计四 用MATLAB进行数值计算五 几类主要模型的MATLAB命令 一 MatLab简介及基本运算 1 MATLAB简介MATLAB Matrix Laboratory 是美国MathWorks公司自20世纪80年代中期推出的数学软件 优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出 到目前为止 其最高版本7 0版已经推出 随着版本的不断升级 它在数值计算及符号计算功能上得到了进一步完善 MATLAB已经发展成为多学科 多种工作平台的功能强大的大型软件 如今 MATLAB已经成为线性代数 自动控制理论 概率论及数理统计 数字信号处理 时间序列分析 动态系统仿真等高级课程的基本教学工具 MATLAB具有用法简易 可灵活运用 程式结构强又兼具延展性 以下为其几个特色 强大的数值计算和工程运算功能先进的资料视觉化功能高阶但简单的程式环境开放及可延伸的架构丰富的程式工具箱 面向专门领域的工具箱 小波工具箱 神经网络工具箱 信号处理工具箱 图像处理工具箱 模糊逻辑工具箱 优化工具箱 鲁棒控制工具箱等几十个不同应用的工具箱 MATLAB能干什么 MATLAB可以进行 数学计算 算法开发 数据采集建模 仿真 原型数据分析 开发和可视化科学和工程图形应用程序的开发 包括图形用户界面的创建 MATLAB广泛应用于 数值计算 图形处理 符号运算 数学建模 系统辨识 小波分析 实时控制 动态仿真等领域 2 MATLAB界面包括 命令窗口 图形窗口 编辑窗口 帮助窗口 3 常用MATLAB系统命令help在线帮助who显示当前变量whos显示当前变量的详细信息clear清空工作间的变量和函数pack整理工作间的内存load把文件调入变量到工作间save把变量存入文件中echo命令回显what显示指定的matlab文件dir显示目录内容 4 基本数学运算 在MATLAB下进行基本数学运算 只需将运算式直接打在提示号 后面 并按Enter键 MATLAB将计算的结果以ans显示 例 求的算术运算结果 1 用键盘在MATLAB指令窗中输入以下内容 12 2 7 4 3 2 2 在上述表达式输入完成后 按 Enter 键 该就指令被执行 3 在指令执行后 MATLAB指令窗中将显示以下结果 ans 2我们也可给运算式的结果设定一个变量x x 5 2 1 3 0 8 10 2 25x 42变量x的值可以在下个语句中调用 y 2 x 1y 85 变量命名规则 1 变量名的大小写敏感 2 变量的第一个字符必须为英文字母 而且不能超过31个字符 3 变量名可以包含下连字符 数字 但不能为空格符 标点 注 数值型变量和符号型变量在MATLAB中是不相同的 它们之间不能直接进行转化 MATLAB提供了一个将数值型转化成符号型的命令 即sym syms 系统预定义的变量ans预设的计算结果的变量名epsMATLAB定义的正的极小值 2 2204e 16pi内建的 值 3 1415926 inf 值 无限大NaN无法定义一个数目 i或j虚数单位i j clear去除所有定义过的变量名称 MATLAB提供基本的算术运算有 加 减 乘 除 幂次方 范例为 5 3 5 3 5 3 5 3 5 3 MATLAB书写表达式的规则与 手写算式 差不多相同 但要求所有表达式都是以纯文本形式输入 如果一个指令过长可以在结尾加上 代表此行指令与下一行连续 例如 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16ans 744 MATLAB常用数学函数三角函数和双曲函数 指数函数 复数函数 其他函数 5 阵列与矩阵MATLAB的运算事实上是以阵列 array 及矩阵 matrix 方式在做运算 阵列强调元素对元素的运算 而矩阵则采用线性代数的运算方式 宣告一变量为阵列或是矩阵时 须用中括号 将元素置于其中 阵列为一维元素所构成 而矩阵为多维元素所组成 例如 x 12345678 一维1x8阵列 例 简单矩阵的输入步骤 1 在键盘上输入下列内容 以 区隔各列的元素 A 1 2 3 4 5 6 7 8 9 2 按 Enter 键 指令被执行 3 在指令执行后 MATLAB指令窗中将显示以下结果 A 123456789 阵列的运算符号除了加减符号外其余的阵列运算符号均须多加 符号 阵列运算功能 注意 一定要多加 符号 加 减 乘 左除 右除 次方 转置矩阵的几种基本变换操作 1 通过在矩阵变量后加 的方法来表示转置运算 a 10 2 12 34 2 4 98 34 6 a ans 10349822341246 2 矩阵求逆 inv a ans 0 01160 0372 0 00150 0176 0 10470 03450 0901 0 0135 0 0045 3 矩阵的特征值 u v eig a u 0 29600 3635 0 3600 0 2925 0 41280 7886 0 9093 0 83520 4985v 48 8395000 19 8451000 10 9943 4 取出上三角和下三角 triu a ans 10212024006 tril a ans 1000342098346 l u lu a l 0 10200 15001 00000 34691 000001 000000u 98 000034 00006 00000 9 79591 91840011 1000 二 MatLab图形绘制功能 基本绘图函数 绘制平面曲线 exp m clf t 0 0 1 3 pi alpha 0 0 1 3 pi plot t sin t r holdon plot alpha 3 exp 0 5 alpha k set gca fontsize 15 fontname timesNewRoman xlabel it t deg ylabel it magnitude title it sinewaveand it Ae alpha itt wave text 6 sin 6 fontsize 15 TheValue it sin t at itt 6 rightarrow bullet HorizontalAlignment right text 2 3 exp 0 5 2 fontsize 15 bullet leftarrowTheValueof it 3e 0 5 it t num2str 3 exp 0 5 2 at it t 2 legend itsin t itAe alphat 注1 num2str string1 num2str string2 用方括号注2 legend请结合图形观察此命令的使用 续 三 MatLab程序设计 1 M文件我们前面所介绍在MATLAB所做的运算 是适合于所要计算的算式不太长或是想以交谈式方式做运算 如果要计算的算式很长有数十行或是须要一再执行的算式 则那样的方式就行不通了 MATLAB提供了所谓的M file的方式 可让使用者自行将指令及算式写成巨集程式然后储存成一个特别的文档 其扩展名是m 譬如picture m 其中的picture就是文件名称 1 编写M脚本文件的步骤点击MATLAB指令窗工具条上的NewFile图标 就可打开如图所示的MATLAB文件编辑调试器 用户即可在空白窗口中编写程序 点击编辑调试器工具条图标 在弹出的Windows标准风格的 保存为 对话框中 选择保存文件夹 键入新编文件名 如area 点击 保存 键 就完成了文件保存 2 运行文件使area m所在目录成为当前目录 系统默认路径 或让该目录处在MATLAB的搜索路径上 然后在指令窗口运行以下指令area例 新建M file area m 计算一个球的体积r input Typeradius输入半径 Area pi r 2 volume 4 3 pi r 3 fprintf 半径Theradiusis 12 5f n r fprintf 面积Theareaofacircleis 12 5f n Area fprintf 体积Thevolumeofasphereis 12 5f n volume 在指令窗口运行指令area 2 M文件函数一个函数M文件与脚本文件类似之处在于它们都是一个有 m扩展名的文本文件 如同脚本M文件一样 函数M文件不进入命令窗口 而是由文本编辑器所创建的外部文本文件 一个函数的M文件与脚本文件在通信方面是不同的 函数与MATLAB工作空间之间的通信 只通过传递给它的变量和通过它所创建的输出变量 在函数内中间变量不出现在MATLAB工作空间 或与MATLAB工作空间不交互 一个函数的M文件的第一行把M文件定义为一个函数 并指定它的名字 它与文件名相同 但没有 m扩展名 它也定义了它的输入和输出变量 M文件函数之间可以互相调用 M文件函数必须遵循以下特定的规则 函数名和文件名必须相同 例如 函数fli存储在名为fli m文件中 例 计算f 1 f 2 f2 3 首先建立一个M函数 fun1 mfunctionY fun1 x Y x 3 2 x 2 x 6 3 x 2 0 05 x 3 14 在指令窗口运行以下指令 fun1 1 fun1 2 fun1 3 fun1 3 ans 12 6023 例 计算阶层函数n fact mfunctionoutput fact n ifn 1output 1 return endoutput n fact n 1 3 MATLAB控制流 for循环结构For循环允许一组命令以固定的和预定的次数重复 For循环的一般形式是 forx array commands end在for和end语句之间的 commands 按数组中的每一列执行一次 for循环可按需要嵌套 例 forn 1 10 x n sin n pi 10 end xx Columns1through70 30900 58780 80900 95111 00000 95110 8090Columns8through100 58780 30900 0000 while循环结构与for循环以固定次数求一组命令的值相反 while循环以不定的次数求一组语句的值 while循环的一般形式是 whileexpression commands end只要在表达式里的所有元素为真 就执行while和end语句之间的 commands 可以利用break命令跳出while循环 while循环可按需要嵌套 例 Fibonacci数组的元素满足Fibonacci规则 且 现要求该数组中第一个大于10000的元素 a 1 1 a 2 1 i 2 whilea i 10000a i 1 a i 1 a i 当现有的元素仍小于10000时 求解下一个元素 i i 1 end i a i i 21ans 10946 if else end分支结构很多情况下 命令的序列必须根据关系的检验有条件地执行 在编程语言里 这种逻辑由某种if else end结构来提供 最简单的if else end结构是 ifexpression commands end如果在表达式中的所有元素为真 非零 那么就执行if和end语言之间的 commands 假如有两个选择 if else end结构是 ifexpressioncommandsevaluatedifTrueelsecommandsevaluatedifFalseend在这里 如果表达式为真 则执行第一组命令 如果表达式是假 则执行第二组命令 当有三个或更多的选择时 If Else End结构采用形式ifexpression1commandsevaluatedifexpression1isTrueelseifexpression2commandsevaluatedifexpression2isTrueelseifexpression3commandsevaluatedifexpression3isTrueelseifexpression4commandsevaluatedifexpression4isTrueelseif elsecommandsevaluatedifnootherexpressionisTrueend最后的这种形式 只和所碰到的 与第一个真值表达式相关的命令被执行 接下来的关系表达式不检验 跳过其余的If Else End结构 而且 最后的else命令可有可无 例 画三次样条函数的图形x 0 0 01 5 fori 1 length x ifx i 0y i x i 3 6 elseifx i 1y i x i 2 3 2 x i 2 2 2 3 elseifx i 2y i x i 2 3 2 x i 2 2 2 3 elseifx i 3y i 4 x i 3 6 elsey i 0 endendplot x y switch case结构例 学生的成绩管理 用来演示switch结构的应用 划分区域 满分 100 优秀 90 99 良好 80 89 及格 60 79 不及格 60 N input 输入分数 switchNcase100 得分为100时S 满分 列为 满分 等级case90 得分在90和99之间S 优秀 列为 优秀 等级case80 得分在80和89之间S 良好 列为 良好 等级case60 得分在60和79之间S 及格 列为 及格 等级otherwise 得分低于60 S 不及格 列为 不及格 等级enddisp S 四 用MatLab进行数值计算 方程求解求解代数方程 组 求解线性方程 求解微分方程 组 sovlefsolvedsovle微分和积分diffint数值积分梯形法trapz二次函数法quadquads 线性优化用命令x lp C A b vlb vub 非线性优化用命令x constr f x0 曲线拟合与插值 五 几类主要模型的MATLAB命令 一 微分方程数值解 1 求微分方程 组 的解析解命令 dsolve 方程1 方程2 方程n 初始条件 自变量 记号 在表达微分方程时 用字母D表示求微分 D2 D3等表示求高阶微分 任何D后所跟的字母为因变量 自变量可以指定或由系统规则选定为确省 2 用Matlab软件求常微分方程的数值解 t x solver f ts x0 options 1 在解n个未知函数的方程组时 x0和x均为n维向量 m 文件中的待解方程组应以x的分量形式写成 2 使用Matlab软件求数值解时 高阶微分方程必须等价地变换成一阶微分方程组 注意 二 解优化问题 1 线性规划有约束极小问题模型 用命令 x fval linprog f A b A1 b1 lb ub 2 非线性规划有约束极小问题模型1 用命令x constr f x0 Examples Findvaluesofxthatminimizef x x1x2x3 startingatthepointx 10 10 10 andsubjecttotheconstraints0 x1 2x2 2x3 72 x1 2x2 2x3 0 x1 2x2 2x3 72 第一步 编写M文件function f g myfun x f x 1 x 2 x 3 g 1 x 1 2 x 2 2 x 3 g 2 x 1 2 x 2 2 x 3 72 第二步 求解在MATLAB工作窗中键入x0 10 10 10 x constr myfun x0 即可 模型2 MATLAB求解此问题的命令是 x fval exitflag output lambda grad hessian fmincon fun x0 A b A1 b1 LB UB nonlcon options p1 p2 fun是目标函数的m 文件名 nonlcon是约束函数C x 和C1 x 的m 文件名 文件输出为 C C1 3 非线性无约束极小问题 用命令x fmin f x0 或用命令x fminu f x0 或用命令x fmins f x0 4 非线性最小二乘问题 用命令x leastsq f x0 或用命令x curvefit f x0 5 二次规划 用命令x qp H c A b 合金强度y与其中含碳量x有密切关系 如下表 根据此表建立y x 并对结果作可信度进行检验 判断x对y影响是否显著 检查数据中有无异常点 由x的取值对y作出预测 2 将17至19岁的运动员每两岁一组分为7组 每组两人测量其旋转定向能力 以考察年龄 x 对这种运动能力 y 的影响 现得到一组数据如下表 合金强度y与其中含碳量x有密切关系 如下表 根据此表建立y x 并对结果作可信度进行检验 判断x对y影响是否显著 检查数据中有无异常点 由x的取值对y作出预测 解 在x y平面上画散点图 直观地知道y与x大致为线性关系 用命令polyfit x y 1 可得y 140 6194x 27 0269 x 0 1 0 01 0 18 x x 0 2 0 21 0 23 y 42 41 5 45 45 5 45 47 5 49 55 50 55 55 5 60 5 polyfit x y 1 plot x y r 三 回归分析 作回归分析用命令 b bint r rint ststs regress y x alpha 可用help查阅此命令的具体用法残差及置信区间可以用rcoplot r rint 画图 x 0 1 0 01 0 18 x x 0 2 0 21 0 23 y 42 41 5 45 45 5 45 47 5 49 55 50 55 55 5 60 5 X ones 12 1 x b bint r rint stats regress y X 0 05 b bint stats rcoplot r rint 结果含义为 0 27 0269 1 140 6194 0的置信区间是 22 3226 31 7313 1的置信区间是 111 7842 169 4546 tji01 m R2 0 9219F 118 0670 p 10 4 R是衡量y与x的相关程度的指标 称为相关系数 R越大 x与y关系越密切 通常R大于0 9才认为相关关系成立 F是一统计指标 p是与F对应的概率 当p 0 05时 回归模型成立 此例中p 0 10 4 0 05 所以 所得回归模型成立 观察所得残差分布图 看到第8个数据的残差置信区间不含零点 此点视为异常点 剔除后重新计算 此时键入 X 8 y 8 b bint r rint stats regress y X b bint stats rcoplot r rint 可以看到 置信区间缩小 R2 F变大 所以应采用修改后的结果 2 将17至19岁的运动员每两岁一组分为7组 每组两人测量其旋转定向能力 以考察年龄 x 对这种运动能力 y 的影响 现得到一组数据如下表 试建立关系y x 并作必要的统计分析 解 在x y平面上画散点图 直观地知道y与x大致为二次函数关系 设模型为y a1x2 a2x a3 此问题可以利用命令polyfit x y 2 来解 也可以象上题一样求解 x 17 2 29 x x x y 20 48 25 13 26 1530 26 1 20 3 19 35 24 35 28 11 26 3 31 4 26 92 25 7 21 3 polyfit x y 2 plot x y r 下面介绍用命令polytool来解 x 17 2 29 x x x y 20 48 25 13 26 15 30 26 1 20 3 19 35 24 35 28 11 26 3 31 4 26 92 25 7 21 3 polytool x y 2 得到一个交互式窗口 首先在命令窗口键入 窗口中绿线为拟合曲线 红线为y的置信区间 可通过移动鼠标的十字线或通过在窗口下方输入来设定x值 窗口左边则输出与x对应的y值及y的置信区间 通过左下方的Export下拉菜单可输出回归系数等 更详细的解释可通过help查阅 3 某厂生产的某产品的销售量与竞争对手的价格x1和本厂的价格x2有关 下表是该产品在10个城市的销售记录 试建立关系y x1 x2 对结果进行检验 若某城市本厂产品售价160 元 对手售价170 元 预测此产品在该城市的销售量 解 这是一个多元回归问题 若设回归模型是线性的 即设y 0 1x1 2x2 那么依然用regress y x alpha 求回归系数 键入 x1 120 140 190 130 155 175 125 145 180 150 x2 100 110 90 150 210 150 250 270 300 250 y 102 100 120 77 46 93 26 69 65 85 x ones 10 1 x1 x2 b bint r rint stats regress y x b bint stats 得b 66 51760 4139 0 2698bint 32 5060165 5411 0 20181 0296 0 4611 0 0785stats 0 65276 57860 0247 p 0 0247 若显著水平取0 01 则模型不能用 R2 0 6527较小 0 1的置信区间包含零点 因此结果不理想 于是设模型为二次函数 此题设模型为纯二次函数 MATLAB提供的多元二项式回归命令为rstool x y model alpha 其中alpha为显著水平 model在下列模型中选一个 linear 线性 purequadratic 纯二次 interaction 交叉 quadratic 完全二次 对此例 在命令窗中键入x 1 rstool x y purequadratic 得到一个对话窗 其意义与前面的对话窗意义类似 若要回答 本厂售价160 对手售价170 预测该市销售量 的问题 只需在下方窗口中分别肩入160和170 就可在左方窗口中读到答案及其置信区间 下拉菜单Export向工作窗输出数据具体操作为 弹出菜单 选all 点击确定 此时可到工作窗中读取数据 可读数据包括 beta 回归系数 rmse 剩余标准差 residuals 残差 本题只要键入 beta rmse residuals 得beta 312 58717 2701 1 7337 0 02280 0037rmse 16 6436residuals 6 6846 12 6703 0 20136 4855 19 65337 9989 11 47375 4303 4 993222 3926大家不妨用此命令选其它模型作一个比较 MATLAB中还包括神经网络工具箱 小波分析工具箱 在网上还可以下载遗传算法工具箱 有兴趣的同学可以借这次机会 结合学习MATLAB 好好学习一下相关理论知识 谢谢大家 最后 祝大家学习 竞赛都取得成功 参考书 数学实验 高等教育出版社参考书 科学计算技术与MATLAB 科学出版社 Thankyouverymuch