I 目录第一章 大学生数学建模竞赛综述 1 第一节 数学建模竞赛的兴起和发展.1 第二节 数学建模竞赛的组织、培训和论文撰写 3 第三节 数学模型的建立过程.4 第四节 几个建模示例.8 第五节 数学建模竞赛题评析.13 第二章 规划与优化模型 20 第一节 奶制品的生产与销售.20 第二节 自来水输送与货机装运.23 第三节 汽车生产与原油采购.27 第四节 接力队选拔和选课策略.31 第五节 饮料厂的生产与检修.35 第三章 微分方程模型 38 第一节 根据规律建模.38 第二节 微元法建模.39 第三节 模拟近似法建模.39 第四节 微分方程模型实例.41 第四章 图论模型 49 第一节 图论的基本概念.49 第二节 图的存储结构.53 第三节 树与图的生成树.54 第四节 网络流及其应用.56 第五节 图的遍历.63 第五章 数学软件使用初步 66 第一节 MATLAB 66 §1.1 MATLAB 简介.66 §1.2 基本运算与函数.66 §1.3 基本平面绘图命令.71 §1.4 基本三维绘图命令.73 §1.5 流程控制.75 §1.6 M 文件 80 §1.7 曲线拟合与插值.81 §1.8 求解非线性规划.87 第二节 Maple.87 §1.1 初识 Maple.88 §1.2 初等数学.89 §1.3 方程求解.92 §1.4 线性代数(矩阵运算)94 §1.5 微积分.97 §1.6 作图.100 II §1.7 函数与过程.103 §1.8 流程控制.107 第三节 LINDO 与LINGO 108 §1.1 LINDO 与LINGO 简介 108 §1.2 用LINDO 求解线性规划(LP)问题 109 §1.3 用LINDO 求解整数规划和二次规划问题.112 §1.4 LINGO 功能简介.114 1 第一章 大学生数学建模竞赛综述 第一节 数学建模竞赛的兴起和发展 一. 数学建模竞赛的历程 ? 20 世纪 60~70 年代进入西方国家的大学(数学建模教材较集中地出现在 70 年代). ? 20 世纪 80 年代初开始进入我国大学; 1987 年出版第 1 本教材(《数学模型》, 姜启源编, 高教 社); 80 年代末估计 30~40 所学校开课(数学系, 讲座). ? 1985 年美国大学生数学建模竞赛开始举办, 1989 年我国大学生开始参加这项竞赛. ? 1992 年我国大学生数学建模竞赛开始举办, 2001 年有 26 省(市、自治区)500 多所学校参加. ? 数学建模竞赛与教学相互促进. 数学模型课已作为数学与应用数学、信息与计算科学等专业 的基础课 二. 广西及我院参赛情况 ? 广西 1994 年开始组织参赛, 当时只有广西大学、桂林电子工业学院参加. ? 2002 年广西已有 13 所本专科院校(1999 年设专科组)53 个队参赛. 但迄今为止参加美国大学生 数模竞赛还是空白. ? 我院 1996 年第一次组队参加全国大学生数学建模竞赛. 获一个三等奖. ? 1998 年再次组队(3 队)参赛, 获广西一等奖一个. ? 1999 年3个队, 全国二等奖 1 个, 广西三等奖 1 个; 2000 年5个队, 全国二等奖 1 个, 广西二、 三等奖各 1 个; 2001 年6个队, 广西三等奖 2 个; 2002 年6个队, 全国二等奖 1 个, 广西一、 二、三等奖各 1 个; 2003 年6个队, 全国二等奖 1 个, 广西二等奖 2 个、三等奖各 1 个. 1. 顺应知识经济时代的潮流, 适应高科技发展的需要教育必须反映并满足社会发展的需求 ? 计算机技术和数学软件的迅速发展, 为数学建模的应用提供了强有力的工具; ? 数学迅速进入一些诸如经济、生态、人口、交通等领域, 为数学建模开拓了许多新的处女地; ? 计算机技术与数学建模在知识经济中起着如虎添翼的作用; ? "数学是一种关键的,普遍的,可应用的技术" 数学"由研究到工业领域的技术转化, 对加强 经济竞争具有重要意义"; ? "计算和建模重新成为中心课题, 它们是数学科学技术转化的主要途径"; ? 教育必须反映并满足社会发展的需求; 2. 培养学生用数学的能力, 探索数学教学改革的途径 ? 数学教育应该培养学生两种能力: "计算数学"(计算、推导、证明…)和"应用数学"(实际 问题建模及模型结果的分析、检验、应用); ? 传统数学教学体系和内容偏重前者, 忽略后者; ? 数学建模引入教学是不打乱现有体系下的教改实验. 数学建模竞赛的迅速发展, 培养了学生创新精神, 推动了高校的教学改革. 三. 数学建模竞赛内容、形式和评奖标准 竞赛内容: 题目由工程技术、管理科学中的实际问题简化而成, 没有事先设定的标准答案, 但留 有充分余地供参赛者发挥其聪明才智和创造精神. 竞赛形式: 三名大学生组成一队, 可以自由地收集资料、调查研究, 使用计算机、互联网和任何 软件, 在三天时间内分工合作完成一篇论文. 评奖标准: 假设的合理性、建模的创造性、结果的正确性和文字表述的清晰程度. 2 大学阶段难得的一次近似于"真刀真枪"的训练, 模拟了毕业后工作时的情况, 既丰富、活跃了 广大同学的课外生活, 也为优秀学生脱颖而出创造了条件. 四. 数学建模竞赛培养学生创新精神, 提高学生综合素质 ? 运用学过的数学知识和计算机(包括选择合适的数学软件)分析和解决实际问题的能力 ? 面对复杂事物的想象力、洞察力、创造力和独立进行研究的能力 ? 关心、投身国家经济建设的意识和理论联系实际的学风 ? 团结合作精神和进行协调的组织能力 ? 勇于参与的竞争意识和不怕困难、奋力攻关的顽强意志 ? 查阅文献、收集资料以及撰写科技论文的文字表达能力 竞赛宗旨: 创新意识 团队精神 重在参与 公平竞争 ? 答卷按省(市、自治区)和全国两级评奖 ? 每年赛题、优秀答卷及获奖名单刊登于次年"数学的实践与认识"和"工程数学"第1期?1999 年的竞赛命名为"99 创维杯全国大学生数学建模竞赛" ? 2000 年的竞赛命名为"2000 网易杯全国大学生数学建模竞赛" ? 2002 年起竞赛命名为"2002 高教社杯全国大学生数学建模竞赛" ? 全国组委会网址: http://csiam.edu.cn/mcm/ ? 竞赛的社会影响不断扩大 五. 数学建模教学和竞赛的教材、辅导材料及参考书 ? 姜启源, 数学模型(第二版), 高等教育出版社, 1993 ? 任善强等, 数学模型(第二版), 重庆大学出版社, 1998 ? 谭永基等, 数学模型, 复旦大学出版社, 1997 ? 刘来福等, 数学模型与数学建模, 北京师范大学出版社, 1997 ? 陈义华, 数学模型, 重庆大学出版社, 1995 ? 徐全智等, 数学建模入门, 电子科技大学出版社, 1996 ? 吴翊等, 数学建模的理论与实践, 国防科技大学出版社, 1999 ? 周义仓等, 数学建模实验, 西安交通大出版社, 1999 ? 白其峥主编, 数学建模案例分析, 海洋出版社, 1999 ? W. F. Lucas 主编, 朱煜民等译, Modules in Applied Mathematics (四卷本, 包括微分方程模型, 政治及有关模型, 离散与系统模型, 生命科学模型), 国防科技大学出版社, 1997(第一卷 1988) ? 姜启源等, 数学实验, 高等教育出版社, 1999 ? 李尚志等, 数学实验, 高等教育出版社, 1999 ? 乐经良等, 数学实验, 高等教育出版社, 1999 ? 谢云逊等, 数学实验, 科学出版社, 1999 ? 叶其孝主编, 大学生数学建模竞赛辅导教材, 湖南教育出版社, 1993 ? 叶其孝主编, 数学建模教育与国际数学建模竞赛, 《工科数学》杂志社, 1994 ? 叶其孝主编, 大学生数学建模竞赛辅导教材(二), 湖南教育出版社, 1997 ? 叶其孝主编, 大学生数学建模竞赛辅导教材(三), 湖南教育出版社, 1998 ? 李大潜主编, 中国大学生数学建模竞赛, 高等教育出版社, 1998 ? 李尚志主编, 数学建模竞赛教程, 江苏教育出版社, 1996 六. 数学建模相关期刊 ? 数学的实践与认识(中, 季刊) ? 工程数学(季刊) 3 ? Mathematical and Computer Modelling (美, 半月刊) ? Applied Mathematical Modelling (美, 月刊) ? The Journal of Undergraduate Mathematics and its Applications(UMAP, 美, 季刊) ? Mathematical Models and Methods in Applied Science(新加坡, 双月刊) 七. 全国数学建模教学与应用会议 第1届(1986, 西安); 第2届(1988, 衡阳); 第3届(1991, 长沙); 第4届(1993, 太原); 第5届(1996, 长春); 第6届(1998, 南昌); 第7届(2000, 郑州); 第8, 9 届(8 月, 北京); 第10 届(2003, 大连) 八. International Conference on the Teaching of Mathematical Modeling and Applications ( ICTMA) ICTMA1(1983)每两年一次, 欧洲; ICTMA8(1997, 澳); ICTMA9(1999, 葡); ICTMA10(2001, 中国北京) 九. 数学建模竞赛中常用的模型 ? 初等数学方法建模(代数、几何、初等概率方法) ? 量纲分析法建模 ? 微分法建模(静态优化模型) ? 微分方程模型(动态模型, 常微部分) ? 差分方程模型 ? 层次分析法建模 ? 随机模型(概率分布方法建模) ? 微分方程模型(偏微部分) ? 稳态模型(稳定性方法建模) ? 图的方法建模(简单的图论方法的应用) ? 逻辑方法建模(合作对策模型等) ? 马氏链模型 ? 随机服务模型 ? 数学规划模型 ? 回归模型 第二节 数学建模竞赛的组织、培训和论文撰写 一. 数学建模竞赛的培训内容 1)建模的基本概念和方法(数学建模课程的主要内容) 2)建模过程中常用的数学方法(微积分、代数、概率外), 主要有: 计算方法(如数值微分和积分、 微分方程数值解、代数方程组解法), 优化方法(如线性、非线性规划), 数理统计(如假设检验、回归 4 分析), 图论(如最短路)等. 只要求知道实际问题与这些数学知识之间的对应关系(如哪些问题可用线性规划求解, 或线性规 划可解决哪些问题), 以及用它们建立模型的方法, 基本上不必涉及模型的求解. 3)合适的数学软件的基本用法. 基本上能完成上述方法的软件, 如Mathematica, Maple 等. 4)历届赛题的研讨. 5)撰写数学建模论文的练习. 其中 1)2)以教师讲授为主, 3)5)以学生实习为主, 4)以学生讨论、教师辅导为主. 二. 数学建模竞赛组队的方式 ? 尽可能地让不同专业的学生组成一队, 以利学科交叉; ? 尽可能地让能力、素质方面不同的学生(创新能力强的, 认真踏实的, 有组织能力的, 文笔好 的, …)组成一队, 以利优势互补; ? 尽可能地让学生在队内充分磨合, 达成默契, 形成"领袖". 三. 数学建模竞赛期间的注意事项 ? 吃透题意, 确定题目; ? 查阅资料、实际调查要适度; ? 把握好用现成的模型和方法, 与自己创新的模型和方法之间的关系; ? 保证基本模型和求解的完成, 在此基础上完善或改进; ? 根据建模的要求, 可以增加、删除甚至修改题目的条件; ? 论文主体由一人完成, 并早些开始写作. 四. 写好论文(答卷)的注意事项 ? 完整——摘要; 问题提出(用自己的语言); 问题分析; 模型假设; 模型建立; 模型求解(算法设 计和计算机实现); 结果(数据、图形); 结果分析和检验(如误差分析、统计检验、灵敏性检验); 优缺 点, 改进方向等, 附录(程序、更多的计算结果、复杂的推导、证明等); ? 摘要——主要模型(名称)、方法和结果, 解决了什么问题, 有何特色等; ? 表述清晰、简明, 给出数学符号的确切含义、模型假设的理由等. 第三节 数学模型的建立过程 一. 什么是数学模型? 我们常见的模型 玩具、照片… —— 实物模型 风洞中的飞机… —— 物理模型 地图、电路图… —— 符号模型 模型是为了一定目的, 对客观事物的一部分进行简缩、抽象、提炼出来的原型的替代物. 模型集中反映了原型中人们需要的那一部分特征. 二. 你碰到过的数学模型——"航行问题" 甲乙两地相距 750 公里, 船从甲到乙顺水航行需 30 小时, 从乙到甲逆水航行需 50 小时, 问船 的速度是多少. 用x表示船速, y 表示水速, 列出方程: ( ) 30 750 ( ) 50 750 x y x y + * = ? * = 求解得到 x=20, y=5, 答: 船速每小时 20 公里. 5 三. 航行问题建立数学模型的基本步骤 ? 作出简化假设(船速、水速为常数); ? 用符号表示有关量(x, y 表示船速和水速); ? 用物理定律(匀速运动的距离等于速度乘以时间)列出数学式子(二元一次方程); ? 求解得到数学解答(x=20, y=5); ? 回答原问题(船速每小时 20 公里). 四. 数学模型(Mathematical Model)和数学建模 Mathematical Modeling) 数学模型: 对于一个现实对象, 为了一个特定目的, 根据其内在规律, 作出必要的简化假设, 运用适当的数学工具, 得到的一个数学结构. 数学建模: 建立数学模型的全过程(包括建立、求解、分析、检验). 五. 数学建模的重要意义 ? 电子计算机的出现及飞速发展. ? 数学以空前的广度和深度向一切领域渗透. 数学建模作为用数学方法解决实际问题的第一步, 越来越受到人们的重视. 六. 建模示例 例1商人们怎样安全过河 问题(智力游戏) 随从们密约, 在河的任一岸, 一旦随从的人数比商人多, 就杀人越货. 但是乘船渡河的方案由商 人决定.商人们怎样才能安全过河? 问题分析 多步决策过程 决策——每一步(此岸到彼岸或彼岸到此岸)船上的人员 要求——在安全的前提下(两岸的随从数不比商人多), 经有限步使全体人员过河 模型构成 xk~第k次渡河前此岸的商人数 yk~第k次渡河前此岸的随从数 xk, yk=0,1,2,3; k=1,2, … sk=(xk , yk)是过程的状态 S 是允许状态集合, S={(x , y)| x=0, y=0,1,2,3; x=3, y=0,1,2,3; x=y=1,2} uk 是第 k 次渡船上的商人数 vk 是第 k 次渡船上的随从数 uk, vk=0,1,2; k=1,2,… dk=(uk , vk) 是决策, D={(u , v)| u+v=1, 2} 是允许决策集合 sk+1=sk+(-1)k dk 是状态转移律 多步决策问题 求dk∈(k=1,2, n), 使sk∈S 按转移律, 由s1=(3,3)到达 sn+1=(0,0). 6 模型求解 1. 穷举法. 编程上机. 2. 图解法. 状态 s=(x,y) ~ 16 个格点, 允许状态 S ~ 10 个圆点, 允许决策 D ~ 移动 1 或2格, k 奇, 左下移; k 偶, 右上移, d1, …d11 给出安全渡河方案. 考虑 4 名商人各带一随从的情况. 例2如何预报人口的增长 世界人口增长概况 年1625 1830 1930 1960 1974 1987 1999 人口(亿) 5 10 20 30 40 50 60 中国人口增长概况 年1908 1933 1953 1964 1982 1990 1995 人口(亿) 3 4.7 6 7 10.1 11.3 12 研究人口变化规律, 控制人口过快增长. 指数增长模型 1. 常用的计算公式 今年人口 x0, 年增长率 r, k 年后人口 xk=x0(1+r)k 2. 马尔萨斯(1788--1834)提出的指数增长模型(1798) x(t)是t时刻的人口, r 是人口(相对)增长率(常数). 由 x t t x t rx t t 解得 0 ( ) rt x t x e = . 设0,(0) dx rx x x dt = = , 得trextx)()(0=trx)1(0+≈.随着时间增加人口按指数规律无限增长. 指数增长模型的应用及局限性 ? 与19 世纪以前欧洲一些地区人口统计数据吻合 ? 适用于 19 世纪后迁往加拿大的欧洲移民后代 7 ? 可用于短期人口增长预测 ? 不符合 19 世纪后多数地区人口增长规律 ? 不能预测较长期的人口增长过程 19 世纪后人口数据是人口增长率 r 不是常数(逐渐下降) 3.阻滞增长模型 (Logistic 模型) 人口增长到一定数量后, 增长率下降的原因: 资源、环境等因素对人口增长的阻滞作用且阻滞 作用随人口数量增加而变大. 假定: r(x)=r-sx(r,s>0), r 是固有增长率(x 很小时), xm 是人口容量(资源、环境能容纳的最大数量, 所以 0 ) ( = m x r , 进而 m r s x = , 所以有 ) 1 ( ) ( m x x r x r ? = . 因为 rx dt dx = , 所以 ) 1 ( ) ( m x x rx x x r dt dx ? = = , 解得 x t x x x e m m rt ( ) ( ) = + ? ? 1 1 0 x(t)~S 形曲线, x 增加先快后慢. 模型的参数估计 用指数增长模型或阻滞增长模型作人口预报, 必须先估计模型参数 r 或r, xm ? 利用统计数据用最小二乘法作拟合 例: 美国人口数据(单位~百万) 1790 1800 1810 1820 1830 …… 1950 1960 1970 1980 3.9 5.3 7.2 9.6 12.9 …… 150.7 179.3 204.0 226.5 r=0.2072, xm=464 ? 专家估计 模型检验 用模型预报 1990 年美国人口, 与实际数据比较. 由]/)1980 ( 1 )[ 1980 ( ) 1980 ( ) 1980 ( ) 1990 ( m x x rx x x x x ? + = ? + = 得5.250 ) 1990 ( = x , 实际为 251.4 (百万). 模型应用——人口预报 用美国 1790~1990 年人口数据重新估计参数, r=0.2083, xm=457.6, 由此得到 x(2000)=275.0, x(2010)=297.9 数学建模的方法和步骤 基本方法: 8 ?机理分析 根据对客观事物特性的认识, 找出反映内部机理的数量规律. ?测试分析 将研究对象看作"黑箱",通过对量测数据 的统计分析, 找出与数据拟合最好的模型 机理. ?二者结合 机理分析建立模型结构,测试分析确定模型参数. 分析没有统一的方法, 主要通过实例研究(Case Studies)来学习. 以下建模主要指机理分析 数学建模的一般步骤: 七. 怎样学习数学建模 数学建模与其说是一门技术, 不如说是一门艺术 技术大致有章可循 艺术无法归纳成普遍适用的准则 数学建模需要有丰富的想象力、敏锐的洞察力、准确的判断力和大胆的创新意识. 开始学习建模时需要学习、分析、评价、改进别人作过的模型, 然后亲自动手, 认真作几个实际 题目. 第四节 几个建模示例 例1. 录像机计数器的用途 问题: 经试验, 一盘录像带从头走到尾, 时间用了 183 分30 秒, 计数器读数从 0000 变到 6152. 在一次使用中录像带已转过大半, 计数器读数为 4580, 问剩下的一段还能否录下 1 小时的节 目? 要求: 不仅回答问题, 而且建立计数器读数与 录像带转过时间的关系. 思考: 计数器读数是均匀增长的吗? 问题分析 录像机计数器的工作原理 录像带运动→右轮盘半径增大→计数器读数增长变慢→录像带运动速度是常数→右轮转速不是 常数. 观察得出: 计数器读数增长越来越慢! 模型准备 模型假设 模型构成 模型求解 模型分析 模型检验 模型应用 9 模型假设 ? 录像带的运动速度是常数 v ; ? 计数器读数 n 与右轮转数 m 成正比, 记m=kn; ? 录像带厚度(加两圈间空隙)为常数 w; ? 空右轮盘半径记作 r ; ? 时间 t=0 时读数 n=0 . 建模目的 建立时间 t 与读数 n 之间的关系(设V, k , w , r 为已知参数) 模型建立 建立 t 与n的函数关系有多种方法 1. 右轮盘转第 i 圈的半径为 r+wi, m 圈的总长度 等于录像带在时间 t 内移动的长度 vt, 所以 ∑ ∞ = = + 1 ) ( 2 n vt wi r π , m=nk 由此得 n v rk n v wk t π π 2 2 2 + = 2. 考察右轮盘面积的变化, 等于录像带厚度乘以转过的长度, 即22[( ) ] r wkn r wvt π + ? = 由此得 2 2 2 wk rk t n n v v π π = + 3. 考察 t 到t+dt 录像带在右轮盘缠绕的长度, 有()2 r wkn kdn vdt π + = 由此得 2 2 2 wk rk t n n v v π π = + 思考 1. 三种建模方法得到同一结果, 但仔细推算会发现稍有差别, 请解释. 2. 模型中有待定参数 r, w, v, k, 确定参数的一种办法是测量或调查, 试设计测量方法. 参数估计 确定参数的另一种方法——测试分析 将模型改记作 2 , t an bn = + 只需估计 a, b 理论上, 已知 t=183.5, n=6152, 再有一组(t, n)数据即可; 实际上, 由于测试有误差, 最好用足够多的数据作拟合. 现有一批测试数据: t 0 20 40 60 80 100 120 140 160 183.5 n 0000 1153 2045 2800 3466 4068 4621 5135 5619 6152 用最小二乘法可得 a=2.51*10-6 , b=1.44*10-2 . 模型检验 应该另外测试一批数据检验模型: 2 t an bn = + , (a=2.51*10-6 , b=1.44*10-2 ) 模型应用 10 1. 回答提出的问题: 由模型算得 n = 4580 时t=118.5 分, 剩下的录像带能录 183.5-118.5 = 65 分钟的节目. 2. 揭示了"t 与n之间呈二次函数关系"这一普遍规律, 当录像带的状态改变时, 只需重新估 计a, b 即可. 例2. 市场经济中的蛛网模型 在自由贸易市场上你注意过这样的现象吗: 一个时期以来某种消费品如猪肉的上市量远大于需 求, 由于销售不畅导致价格下降, 生产者发现养猪赔钱, 于是转而经营其它农副业. 过一段时间猪肉 上市量就会大减, 供不应求将导致价格上涨. 生产者看到有利可图, 又重操旧业, 这样下一个时期会 重现供大于求、价格下降的局面. 在没有外界干预的情况下, 这种现象将如此循环下去. 在完全自由竞争的市场经济中上述现象通常是不可避免的. 因为商品的价格呈由消费者的需求 关系决定的, 商品数量越多价格越低. 而下一时期商品的数量由生产者的供应关系决定, 商品价格 越低生产的数量就越少. 这样的需求和供应关系决定了市场经济中商品的价格和数量必然是振荡的. 在现实世界里这样的振荡会出现不同的形式, 有的振幅渐小趋向平稳, 有的则振幅越来越大, 如果 没有外界如政府的干预, 将导致经济崩溃. 请你对上述现象进行分析, 给出市场经济趋于稳定的条件. 并建立数学模型, 对结果进行解释, 并讨论当市场经济不稳定时政府可以采取什么样的干预措施, 并对你的模型作适当推广. 问题分析 现象 问题 ? 描述商品数量与价格的变化规律. ? 商品数量与价格的振荡在什么条件下趋向稳定. ? 当不稳定时政府能采取什么干预手段使之稳定. 模型建立 xk——第k时段商品数量; yk——第k时段商品价格 消费者的需求关系→需求函数 yk=f(xk)减函数. 生产者的供应关系→供应函数 xk+1=h(yk), yk=g(xk+1)增函数. f 与g的交点 P0(x0,y0)——平衡点,一旦 xk=x0,则yk=y0, xk+1,xk+2,…=x0, yk+1,yk+2,…=y0. 设x1 偏离 x0, x1→y1→x2→y2→x3→…, xk→x0, yk→y0, P1→P2→P3→…→P0. 11 KfKg, P0 是不稳定平衡点 在P0 点附近用直线近似曲线 ( ) k k y f x = ? 0 0 ( ) ( 0) k k y y x x α α 1 ( ) k k x h y + = ? 1 0 0 ( ) ( 0) k k x x y y β β 由)(001xxxxkk??=?+αβ 得)()(0101xxxxkk??=?+αβ 当1<αβ 时, 1 ( ) ( ) f g K K α β = < = , 0 x xk → , P0 稳定 当1>αβ 时, 1 ( ) ( ) f g K K α β → k x , P0 不稳定 结果解释 考察 α, β 的含义 0 0 ( ), k k y y x x α ? = ? ? α——商品数量减少 1 单位, 价格上涨幅度 1 0 0 ( ), k k x x y y β + ? = ? β——价格上涨 1 单位, (下时段)供应的增量 α——消费者对需求的敏感程度, α 小, 有利于经济稳定 β——生产者对价格的敏感程度, β 小, 有利于经济稳定 由此 αβ<1, 经济稳定 经济不稳定时政府的干预办法 1. 使α尽量小, 如α=0, 需求曲线变为水平, 以行政手段控制价格不变 2. 使β尽量小, 如β=0, 供应曲线变为竖直, 靠经济实力控制数量不变 例3. 报童模型 报童每天清晨从报站批发报纸零售, 晚上将没有卖完的报纸退回. 设每份报纸的批发价为 b, 零 售价为 a, 退回价为 c, 且设 a >b > c>0. 因此, 报童每出售一份报纸赚钱(a-b), 退回一份报纸赔(b-c), 报童每天如果批发太少, 不够卖的话就会少赚钱; 如果批发的报纸太多, 卖不完的话就会赔钱. 报童 12 应如何确定他每天批发的报纸的数量, 才能获得最大的收益? 问题分析 购进太多→卖不完退回→赔钱 购进太少→不够销售→赚钱少 应根据需求确定购进量, 是否存在一个合适的购进量? 每天需求量是随机的, 所以每天收入是随机的. 优化问题的目标函数应是长期的日平均收入, 等于每天收入的期望. 调查需求量的随机规律——每天需求量为 r 的概率 f(r), r=0,1,2…. 模型建立 ? 设每天购进 n 份, 日平均收入为 G(n). ? 已知售出一份赚 a-b; 退回一份赔 b-c. r≤n→售出 r, 退回 n-r, 赚(a-b)r, 赔(b-c)(n-r) r>n→售出 n, 赚(a-b)n 0 1 n r r n G n a b r b c n r f r a b nf r ∞ = = + ∑ ∑ 求n使G(n)最大. 模型求解 将r视为连续变量, f(r)=p(r)(概率密度), 0 n n G n a b r b c n r p r dr a b np r dr ∞ ∫ ∫ . 由上可得 0 n dG a b np n b c p r dr a b np n dn n a b p r dr ∞ + ? ∫ 0 n n b c p r dr a b p r dr ∞ ∫ ∫ . 令0dG dn = 得0()()nnprdr a b b c p r dr ∞ ? = ? ∫ ∫ . 结果解释 1 2 0 n n p r dr P p r dr P ∞ = = ∫ ∫ , 取n使12PabPbc?=?其中 a-b ~售出一份赚的钱, b-c ~退回一份赔的钱. a b n b c n 13 第五节 数学建模竞赛题评析 例1投资的收益和风险(1998 年全国大学生数学建模竞赛题目 A 题) 市场上有 n 种资产(如股票、债券、…)Si(i=1, …, n)供投资者选择, 某公司有数额为 M 的一 笔相当大的资金可用作一个时期的投资. 公司财务分析人员对这 n 种资产进行了评估, 估算出在这 一时期内购买 Si 的平均收益率为 ri, 并预测出购买 Si 的风险损失为 qi. 考虑到投资越分散, 总的风 险越小, 公司确定, 当用这笔资金购买若干种资产时, 总体风险可用所投资的 Si 中最大的一个风险 来度量. 购买 Si 要付交易费, 费率为 pi, 并且当购买额不超过给定值 ui 时, 交易费按购买 ui 计算(不买当 然无须付费). 另外, 假定同期银行存款利率是 r0, 且既无交易费又无风险. (r0=5%) 1) 以知 n=4 时的相关数据如下: Si ri(%) qi(%) pi(%) ui(元) S1 28 2.5 1 103 S2 21 1.5 2 198 S3 23 5.5 4.5 52 S4 23 2.6 6.5 40 试给该公司设计一种投资组合方案, 即用给定的资金 M, 有选择地购买若干种资产或存银行生息, 使净收益尽可能大, 而总体风险尽可能小. 2) 试就一般情况对以上问题进行讨论, 并利用以下数据进行计算. Si ri(%) qi(%) pi(%) ui(元) S1 9.6 42 2.1 181 S2 18.5 54 3.2 407 S3 49.4 60 6.0 428 S4 23.9 42 1.5 549 S5 8.1 1.2 7.6 270 S6 14 39 3.4 397 S7 40.7 68 5.6 178 S8 31.2 33.4 3.1 220 S9 33.6 53.3 2.7 475 S10 36.8 40 2.9 248 S11 11.8 31 5.1 195 S12 9 5.5 5.7 320 S13 35 46 2.7 267 S14 9.4 5.3 4.5 328 S15 15 23 7.6 131 问题分析 这是一个优化问题, 可设决策变量为每种资产的投资额(投资组合), 目标是使得净收益最大, 且 整体风险最小. 但是二者矛盾. 可采取下面几种方案考虑: ? 在一定风险下收益最大的决策. ? 在一定收益下风险最小的决策. ? 收益和风险按一定比例组合最优的决策. 14 这时可得到多组解, 由决策者进行选择. 冒险型投资者从中选择高风险下收益最大的决策, 而 保守型投资者则可从低风险下的决策中选取. 模型建立 用数学符号和式子表述决策变量、构造目标函数、确定约束条件. 决策变量: xi—对Si(i=0,1,…n)的投资, S0 表示存入银行. 目标函数: ? 总收益—投资 Si 的净收益减去交易费, 对i求和. ? 总体风险—投资 Si 的风险, 对i求最大值. 约束条件: 对Si 的投资 xi 加交易费, 对i求和, 不超过给定资金 M. 1) 投资 Si 的交易费、净收益、风险、资金表达式 交易费 0, 0 ( ) , 0 , i i i i i i i i i i i x c x p u x u p x x u = ? ? = < < ? ? ≥ ? 图形如下: 净收益(收益率 ri) ( ) ( ) i i i i i i R x r x c x = ? 风险(风险损失率 qi) ( ) i i i i Q x q x = 资金 0 0 0 0,1,0 i i i i i f x x c x i n c x q 2) 投资方案总收益、总体风险、资金表达式 投资方案 0 1 ( , , ) n x x x x = " 总收益 0 ( ) ( ) n i i i R x R x = = ∑ 总体风险 0 ( ) max ( ) i i i n Q x Q x ≤ ≤ = 资金 15 0 ( ) ( ) n i i i F x f x = = ∑ 3)两目标(总收益、总体风险)优化模型 ( ) min ( ) , 0 ( ) x Q x F x M x R x ? ? ? ? ? ? = ≥ ? ? ? ? ? ? ? ? ? ? ? 4)单目标优化模型 模型 M1: 确定风险水平 q , 记Mqk=,max R(x) s.t. Q(x)≤k F(x)=M, x≥0 模型 M2: 确定盈利水平 r , 记hrM = , min Q(x) s.t. R(x) ≥h F(x)=M, x≥0 模型 M3: 确定投资者对风险—收益的相对偏好参数 ρ>0, min S(x)= ρQ(x)-(1-ρ)R(x) s.t. F(x)=M, x≥0 模型简化 交易费 ci(xi)当ui 很小且 M 很大时, 有ci(xi)=pixi. 对资金约束, 由0()()niiiFxfx= i i i i i f x x c x = + , F(x)=M 得到 0 (1 ) n i i i p x M = + = ∑ 设M=1, 投资 Si 的比例 (1 ) , 0,1,..., i i i y p x i n = + = 由此, 模型 M1 简化为 max 0 ( ) n i i i i r p x = ? ∑ s.t. , 1,2,..., i i q x k i n ≤ = 0 (1 ) 1, 0 n i i i p x x = + = ≥ ∑ 这是一个线性规划模型(LP1). 模型 M2 简化为 min ) ( max 0 i i n i x q ≤ ≤ s.t. h x p r n i i i i ≥ ? ∑ =0 ) ( 16 0 , 1 ) 1 ( 0 ≥ = + ∑ = x x p n i i i 这是一个极大极小规划模型. 引入人工变量 xn+1, 上面模型可化为线性规划模型(LP2): min xn+1 s.t. qixi≤xn+1, i=1,…,n h x p r n i i i i ≥ ? ∑ =0 ) ( 0 , 1 ) 1 ( 0 ≥ = + ∑ = x x p n i i i 模型 M3 的简化为 min ∑ = + ? ? ? = n i i i i n x p r p x x L 0 1 ) ( ) 1 ( ) ( ρ s.t. 0 , 1 ) 1 ( 0 ≥ = + ∑ = x x p n i i i qixi≤xn+1, i=1,…,n 这是一个线性规划模型(LP3). 模型求解 LP1, LP2, LP3 都很容易用 MATLAB, MATHEMATICA 或其它数学软件求解. LP1 求解结果: 风险水平取 k=0~2.5%, 得投资比例 y0~y4. 风险 K(%) 收益 R(%) y0 y1 y2 y3 y4 0 5.0000 1.0000 0 0 0 0 0.1 7.5528 0.8316 0.0404 0.0680 0.0190 0.0410 0.2 10.1055 0.6633 0.0808 0.1360 0.0380 0.0819 0.3 12.6583 0.4949 0.1212 0.2040 0.0570 0.1229 0.4 15.2110 0.3266 0.1616 0.2720 0.0760 0.1638 0.5 17.7638 0.1582 0.2020 0.3400 0.0950 0.2048 0.6 20.1908 0 0.2424 0.4080 0.1140 0.2356 0.7 20.6607 0 0.2828 0.4760 0.1330 0.1082 0.8 21.1243 0 0.3232 0.5440 0.1328 0 0.9 21.5520 0 0.3636 0.6120 0.0244 0 1.0 21.9020 0 0.4040 0.5960 0 0 1.5 23.5392 0 0.6060 0.3940 0 0 2.0 25.1765 0 0.8080 0.1920 0 0 2.1 25.5039 0 0.8484 0.1516 0 0 2.2 25.8314 0 0.8888 0.1112 0 0 2.3 26.1588 0 0.9292 0.0708 0 0 2.4 26.4863 0 0.9696 0.0304 0 0 2.5 26.7327 0 1.0000 0 0 0 LP1 风险 K 和收益 R 的关系 17 对于风险和收益没有特殊偏好的投资者来说, 应该选择图中曲线的拐点(K*, R*), 大约是 K*=0.6%, R*=20%. LP3 求解结果: 偏好系数 ρ 取0.76~0.97. ρ 风险 K(%) 收益 R(%) y0 y1 y2 y3 y4 0.76 2.1752 36.7327 0 1.0000 0 0 0 0.77 0.9225 21.6482 0 0.3727 0.6273 0 0 0.81 0.9225 21.6482 0 0.3727 0.6273 0 0 0.82 0.7849 21.0599 0 0.3171 0.5338 0.1494 0 0.83 0.5940 20.1624 0 0.2400 0.4039 0.1129 0.2432 0.96 0.5940 20.1624 0 0.2400 0.4039 0.1129 0.2432 0.97 0.0000 5.0000 1.0000 0 0 0 0 ρ 与风险 K 的关系 ρ 与收益 R 的关系 风险 K 和收益 R 的关系 结果与 LP1 相同. 例2节水洗衣机(1996 年全国大学生数学建模竞赛题目 B 题) 我国淡水资源有限, 节约用水人人有责. 洗衣在家庭用水中占有相当大的份额, 目前洗衣机已 非常普及, 节约洗衣机用水十分重要. 假设在放人衣物和洗涤剂后洗衣机的运行过程为: 加水—漂洗—脱水—加水—漂洗—脱水—…—加水—漂洗—脱水(称"加水—漂洗—脱水"为运行一轮). 请 为洗衣机设计一种程序(包括运行多少轮、 每轮加水量等), 使得在满足一定洗涤效果的条件下)总用 水量最少. 选用合理的数据进行计算. 对照目前常用的洗衣机的运行情况, 对你的模型和结果作出 评价. 问题分析 ? 不要求对洗衣的微观机制(物理、化学方面)深入研究, 只需从宏观层次去把握. ? 洗衣的基本原理(宏观上)是用洗涤剂通过漂洗把吸附在衣物上的污物溶于水中, 再脱去污水 18 带走污物. ? 洗衣的过程是通过"加水——漂洗——脱水" 程序的反复运行, 使残留在衣物上的污物越来越 少, 直到满意的程度. ? 洗涤剂也是不希望留在衣物上的东西, 可将"污物"定义为衣物上原有污物与洗涤剂的总和. 基本假设 ? 每轮漂洗后污物均匀地溶于水中; ? 每轮脱水后衣物含水量为常数 c; 基本模型 x0——初始污物量 uk——第k轮加水量; xk——第k轮脱水后污物量(k=1,2,…) 每轮脱水前后污物在水中的浓度不变 c x u x 1 1 0 = , c x c u x 2 2 1 = + ,…, c x c u x n n n = + ? ? 1 1 由此得 ) ( ) ( 2 1 0 c u c u u c x x n n n + + = " 在最终污物量与初始污物量之比 xn/x0 小于给定的 ε(清洁度)条件下, 求各轮加水量 uk (k=1,…n), 使总用水量最小. 由此得到模型 1 min k n k u k u = ∑ s.t. 1 2 ( ) ( ) n n c u u c u c ε < + + " . 即12min ( ) ( ) k n u u u c u c " s.t. 1 2 ( ) ( ) n u u c u c a + + = " . 由于几何平均值小于 (等于)算术平均值, 得12nuucuc " . 所以第 2~n 轮加水量 uk= u(常数), 第1轮加水量 u1= u +c. 模型简化为 min nu c + s.t. ( )n c u c ε < + . 又u=cx, 进一步有 , min n u nx s.t. 1 ( ) 1 n x ε < + . 而19 , ln ln(1 ) n x α α ε = = ? + . 得min ln(1 ) x x x α + . 所以 x→0. u=cx→0, n→∞. 考虑加水量限制21vuv≤≤,得min max n n n ≤ ≤ . 其中min 2 [ ] 1 ln(1 / ) n v c α = + + , max 1 [ ] 1 ln(1 / ) n v c α = + + . 作进一步的考虑. m 是衣物重量, 每轮脱水后衣物含水量 c=am, 浸泡衣物所需水量 d=bm, 洗衣 机最小注入水量 v0, 洗衣机水量下限 v1=v0+bm. 用以下实验数据测定 a, b. m(kg) 0.5 1.0 1.5 2.0 2.5 3.0 3.5 c(l) 0.30 0.61 0.92 1.20 1.42 1.88 2.15 d(l) 2.56 5.02 7.48 9.87 12.3 15.7 18.6 得到 a=0.6, b=5.0. 模型结果(单位: 衣物重量(kg); 用水量 (l)) 衣物重量 洗衣次数 第1次用水量 第2次用水量 第3次用水量 总用水量 清洁度 ε 1.0 2 29.0 27.4 56.4 0.003 1.2 2 30.1 28.1 58.2 0.004 1.4 2 31.7 29.4 61.1 0.005 1.6 2 36.2 33.6 69.8 0.005 1.8 2 40.7 37.9 78.6 0.005 2.0 3 34.3 31.1 31.1 96.5 0.0008 2.2 3 35.4 31.8 31.8 99.0 0.0009 2.4 3 36.4 32.6 32.6 101.6 0.0011 2.6 3 37.5 33.3 33.3 104.1 0.0013 2.8 3 38.5 34.1 34.1 106.7 0.0015 3.0 3 39.6 34.8 34.8 109.2 0.0017 与实际比较: 海棠洗衣机 XQB42-1 型, 3kg 衣物用水 127 升. 本模型可节约 18 升. 20 第二章 规划与优化模型 实际问题中的优化模型可归结为: min(or max) z=f(x), x=(x1,…,xn)T s.t. gi(x)≤0,(=0 或≥0) i =1,…,m x 称为决策变量, f(x)称为目标函数, gi(x)≤0 称为约束条件. 具有这样形式的模型称为规划模型. 数学规划可分为: 1. 线性规划——f(x) 、gi(x)均为线性函数, 且xi≥0. 2. 非线性规划——f(x)或某 gi(x)为非线性函数. 3. 整数规划——xi ∈Z . 第一节 奶制品的生产与销售 一奶制品加工厂用牛奶生产 Al, A2 两种奶制品, 1 桶牛奶可以在设备甲上用 12 小时加工成 3 公斤 A1, 或者在设备乙上用 8 小时加工成 4 公斤 A2. 根据市场需求, 生产的 A1, A2 全部能售出, 且每公斤 A1 获利 24 元, 每公斤 A2 获利 16 元. 现在加工厂每天能得到 50 桶牛奶的供应, 每天正式工人总的劳 动时间为 480 小时, 并且设备甲每天至多能加工 100 公斤 A1, 设备乙的加工能力没有限制. 试为该厂 制订一个生产计划, 使每天获利最大, 并进一步讨论以下 3 个问题: 1)若用 35 元可以买到 1 桶牛奶, 应否作这项投资?若投资, 每天最多购买多少桶牛奶? 2)若可以聘用临时工人以增加劳动时间, 付给临时工人的工资最多是每小时几元? 3)由于市场需求变化, 每公斤 A1 的获利增加到 30 元, 应否改变生产计划? 问题分析 空间层次 工厂级: 根据外部需求和内部设备、人力、原料等条件, 以最大利润为目标制订产品生产计划; 车间级: 根据生产计划、工艺流程、资源约束及费用参数等, 以最小成本为目标制订生产批量计 划. 时间层次 若短时间内外部需求和内部资源等不随时间变化, 可制订单阶段生产计划, 否则应制订多阶段 生产计划. 12 小时, 3 公斤 A1 获利 24 元/公斤 1 桶牛奶 或8小时, 4 公斤 A2 获利 16 元/公斤 每天: 50 桶牛奶, 时间 480 小时, 至多加工 100 公斤 A1. 制订生产计划, 使每天获利最大 ? 35 元可买到 1 桶牛奶, 要买吗. 若买, 每天最多买多少? ? 可聘用临时工人, 付出的工资最多是每小时几元? ? A1 的获利增加到 30 元/公斤, 应否改变生产计划? 建立线性规划模型(LP) 决策变量 x1 桶牛奶生产 A1, 获利 24*3x1, x2 桶牛奶生产 A2, 获利 16*4x2 目标函数 每天获利 max z=72x1+64x2 21 约束条件 原料供应 x1+x2≤50 劳动时间 12x1+8x2≤480 加工能力 3x1≤100 非负约束 x1,x2≥0 模型分析与假设 ? 比例性 xi 对目标函数的"贡献"与xi 取值成正比. xi 对约束条件的"贡献"与xi 取值成正比. ? 可加性 xi 对目标函数的"贡献"与xj 取值无关. xi 对约束条件的"贡献"与xj 取值无关. ? 连续性 xi 取值连续 线性规划模型 ? A1, A2 每公斤的获利是与各自产量无关的常数. ? 每桶牛奶加工出 A1, A2 的数量和时间是与各自产量无关的常数. ? A1, A2 每公斤的获利是与相互产量无关的常数. ? 每桶牛奶加工出 A1, A2 的数量和时间是与相互产量无关的常数. ? 加工 A1, A2 的牛奶桶数是实数. 模型求解 约束条件 x1+x2≤50 L1: x1+x2=50 12x1+8x2≤480 L2: 12x1+8x2=480 3x1≤100 L3: 3x1=100 x1,x2≥0 L4:, x1=0, L5 : x2=0 目标函数: max z=72x1+64x2 1. 图解法 z=c (常数) 是等值线, 目标函数和约束条件是线性函数, 可行域为直线段围成的凸多边形, 目标 函数的等值线为直线, 最优解一定在凸多边形的某个顶点取得. 结果: 在B(20,30)点得到最优解. 图解法直观易行, 但当变量和约束条件较多时无法应用. 2. 用LINDO 6.1 软件求解 max 72x1+64x2 st 2) x1+x2<50 22 3) 12x1+8x2<480 4) 3x1<100 end 求LINDO 求解结果(要求作灵敏性分析): LP OPTIMUM FOUND AT STEP 2 OBJECTIVE FUNCTION VALUE 1) 3360.000 VARIABLE VALUE REDUCED COST X1 20.000000 0.000000 X2 30.000000 0.000000 ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 48.000000 3) 0.000000 2.000000 4) 40.000000 0.000000 NO. ITERATIONS= 2 RANGES IN WHICH THE BASIS IS UNCHANGED: OBJ COEFFICIENT RANGES VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE X1 72.000000 24.000000 8.000000 X2 64.000000 8.000000 16.000000 RIGHTHAND SIDE RANGES ROW CURRENT ALLOWABLE ALLOWABLE RHS INCREASE DECREASE 2 50.000000 10.000000 6.666667 3 480.000000 53.333332 80.000000 4 100.000000 INFINITY 40.000000 结果: 20 桶牛奶生产 A1, 30 桶生产 A2, 利润 3360 元. 结果解释 三种资源: 原料无剩余, 时间无剩余, 加工能力剩余 40. "资源"剩余为零的约束为紧约束(有效约束), 最优解下"资源"增加 1 单位时"效益"的增量——影 子价格, 原料增加 1 单位, 利润增长 48 时间增加 1 单位, 利润增长 2, 加工能力增长不影响利润. ? 35 元可买到 1 桶牛奶, 要买吗? 35 <48, 应该买! ? 聘用临时工人付出的工资最多每小时几元? 2 元! 最优解不变时目标函数系数允许变化范围(约束条件不变) x1 系数范围(64,96) x2 系数范围(48,72) x1 系数为 30*3=90, 在允许范围内 ? A1 的获利增加到 30 元/公斤, 应否改变生产计划? 不变! 影子价格有意义时约束右端的允许变化范围 (目标函数不变) 原料最多增加 10 时间最多增加 53 ? 35 元可买到 1 桶牛奶, 每天最多买多少?最多买 10 桶! 23 思考题 前面给出的 A1, A2 两种奶制品的生产条件、利润, 及工厂的"资源"限制全都不变. 为增加工厂 的获利, 开发了奶制品的深加工技术: 用2小时和 3 元加工费, 可将 1 公斤 A1 加工成 0.8 公斤高级奶 制品 B1, 也可将 1 公斤 A2 加工成 0.75 公斤高级奶制品 B2, 每公斤月, 能获利 44 元, 每公斤 B2 能获 利32 元. 试为该厂制订一个生产销售计划, 使每天的净利润最大, 并讨论以下问题: 1)若投资 30 元可以增加供应 1 桶牛奶, 投资 3 元可以增加 1 小时劳动时间, 应否作这些投资?若 每天投资 150 元, 可赚回多少? 2)每公斤高级奶制品 B1, B2 的获利经常有 10%的波动, 对制订的生产销售计划有无影响?若每公 斤B1 的获利下降 10%, 计划应该变化吗? 第二节 自来水输送与货机装运 运输问题一般是指生产、 生活物资从若干供应点运送到一些需求点, 各种类型的货物装箱, 由于 受体积、重量等限制, 如何搭配装载, 使获利最高, 或装箱数量最少, 怎样安排输送方案使运费最小, 或利润最大; 例1自来水输送 某市有甲、乙、丙、丁四个居民区, 自来水由 A, 扫, C 三个水库供应. 四个区每天必须得到保证 的基本生活用水量分别为 30, 70, 10, 10 千吨, 但由于水源紧张, 三个水库每天最多只能分别供应 50, 60, 50 千吨自来水. 由于地理位置的差别, 自来水公司从各水库向各区送水所需付出的引水管理费不 同(见表, 其中 C 水库与丁区之间没有输水管道), 其他管理费用都是 450 元/千吨. 根据公司规定, 各区用户按照统一标准900元/千吨收费. 此外, 四个区都向公司申请了额外用水量, 分别为每天50, 70, 20, 40 千吨. 该公司应如何分配供水量, 才能获利最多? 为了增加供水量, 自来水公司正在考虑进行水库改造, 使三个水库每天的最大供水量都提高一 倍, 问那时供水方案应如何改变?公司利润可增加到多少? 引水管理费(元/千吨) 甲乙丙丁A160 130 220 170 B 140 130 190 150 C 190 200 230 / 问题分析 总供水量: 160< 总需求量: 120+180=300. 收入: 900 元/千吨, 总收入 900*160=144,000(元). 支出: 引水管理费. 其他费用: 450 元/千吨, 其他支出: 450*160=72,000(元). 确定送水方案使利润最大 使引水管理费最小. 模型建立 确定 3 个水库向 4 个小区的供水量 决策变量: 水库 i 向j区的日供水量为 xij(x34=0). 目标函数: min z=160x11+130x12+220x13+170x14+140x21+130x22+190x23+150x24+190x31+200x32+230x33 s.t. x11+x12+x13+x14=50 x21+x22+x23+x24=60 x31+x32+x33 =50 24 30≤x11+x21+x31≤80 70≤x12+x22+x32≤140 10≤x13+x23+x33≤30 10≤x14+x24 ≤50 模型求解 用LINDO 求解以上模型(需要另外编写 LINGO 模型), 得OBJECTIVE FUNCTION VALUE 1) 24400.00 VARIABLE VALUE REDUCED COST X11 0.000000 30.000000 X12 50.000000 0.000000 X13 0.000000 50.000000 X14 0.000000 20.000000 X21 0.000000 10.000000 X22 50.000000 0.000000 X23 0.000000 20.000000 X24 10.000000 0.000000 X31 40.000000 0.000000 X32 0.000000 10.000000 X33 10.000000 0.000000 结果解释 引水管理费: 24400(元) 利润=总收入- 其它费用 - 引水管理费 =144000-72000-24400 = 47600(元) 问题讨论 每个水库最大供水量都提高一倍. 总供水量(320) > 总需求量(300), 确定送水方案使利润最大. 利润= 收入(900) – 其它费用(450) –引水管理费 利润(元/千吨) 甲乙丙丁A290 320 230 280 B 310 320 260 300 C 260 250 220 / 目标函数 max z=290x11+320x12+230x13+280x14+310x21+320x22+260x23+300x24+260x31+250x32+220x33 25 A: x11+x12+x13+x14=50 x11+x12+x13+x14≤100 B, C 类似处理, 需求约束可以不变 求解结果 OBJECTIVE FUNCTION VALUE 1) 88700.00 VARIABLE VALUE REDUCED COST X11 0.000000 20.000000 X12 100.000000 0.000000 X13 0.000000 40.000000 X14 0.000000 20.000000 X21 30.000000 0.000000 X22 40.000000 0.000000 X23 0.000000 10.000000 X24 50.000000 0.000000 X31 50.000000 0.000000 X32 0.000000 20.000000 X33 30.000000 0.000000 结果解释 总利润 88700(元) 例2货机装运 三个货舱最大载重(吨), 最大容积(米3) . 飞机平衡条件 三个货舱中实际载重必须与其最大载重成比例 重量(吨) 空间(米3/吨) 利润(元/吨) 货物 1 18 480 3100 货物 2 15 650 3800 货物 3 23 580 3500 货物 4 12 390 2850 问如何装运, 使本次飞行获利最大? 模型假设 1、每种货物可以分割到任意小; 2、每种货物可以在一个或多个货舱中任意分布; 3、多种货物可以混装, 并保证不留空隙; 4、所给出的数据都是精确的, 没有误差. 模型建立 26 决策变量 xij--第i种货物装入第 j 个货舱的重量(吨) i=1,2,3,4, j=1,2,3 (分别代表前、中、后仓) 目标函数(利润) max Z=3100(x11+x12+x13)+3800(x21+x22+x23) +3500(x31+x32+x33)+2850(x41+x42+x43) 约束条件 xij--第i种货物装入第 j 个货舱的重量 x11+x21+x31+x41≤10 x12+x22+x32+x42≤16 x13+x23+x33+x43≤8 480x11+650x21+580x31+390x41≤6800 480x12+650x22+580x32 +390x42≤8700 480x13+650x23+580x33 +390x43≤5300 (x11+x21+x31+x41)/10= (x12+x22+x32+x42)/16 = (x13+x23+x33+x43)/8 x11+x12+x13≤18 x21+x22+x23≤15 x31+x32+x33 ≤23 x41+x42+x43 ≤12 模型求解 OBJECTIVE FUNCTION VALUE 1) 121515.8 VARIABLE VALUE REDUCED COST X11 0.000000 400.000000 X12 0.000000 57.894737 X13 0.000000 400.000000 X21 10.000000 0.000000 X22 0.000000 239.473679 X23 5.000000 0.000000 X31 0.000000 0.000000 X32 12.947369 0.000000 X33 3.000000 0.000000 X41 0.000000 650.000000 X42 3.052632 0.000000 X43 0.000000 650.000000 货物 2: 前仓 10,后仓 5; 货物 3: 中仓 13, 后仓 3; 货物 4: 中仓 3. 最大利润约 121516 元27 第三节 汽车生产与原油采购 例1汽车厂生产计划 一汽车厂生产小、中、大三种类型的汽车,已知各类型每辆车对钢材、劳动时间的需求,利润 以及每月工厂钢材、劳动时间的现有量如表中所示.试制订月生产计划,使工厂的利润最大. 进一步讨论:由于各种条件限制,如果生产某一类型汽车,则至少要生产 80 辆,那么最优的生 产计划应作何改变. 小型 中型 大型 现有量 钢材(吨) 1.5 3 5 600 劳动时间(小时) 280 250 400 60000 利润(万元) 2 3 4 模型建立 设每月生产小、中、大型汽车的数量分别为 x1, x2, x3 , 则max z=2x1+3x2+4x3 s.t. 1.5x1+3x2+5x3≤600 280x1+250x2+400x3≤60000 x1,x2,x3 ≥0 模型求解 OBJECTIVE FUNCTION VALUE 1) 632.2581 VARIABLE VALUE REDUCED COST X1 64.516129 0.000000 X2 167.741928 0.000000 X3 0.000000 0.946237 ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 0.731183 3) 0.000000 0.003226 结果为小数, 怎么办? 1)舍去小数: 取x1=64, x2=167, 算出目标函数值 z=629, 与LP 最优值 632.2581 相差不大. 2)试探: 如取 x1=65, x2=167; x1=64, x2=168 等, 计算函数值 z, 通过比较可能得到更优的解. 但必 须检验它们是否满足约束条件. 为什么? 3)模型中增加条件: x1, x2, x3 均为整数, 重新求解. 模型求解 整数规划(Integer Programming, 简记 IP) max z=2x1+3x2+4x3 s.t. 1.5x1+3x2+5x3≤600 280x1+250x2+400x3≤60000 x1,x2,x3 为非负整数 IP可用LINDO直接求解: max 2x1+3x2+4x3 st 1.5x1+3x2+5x3<600 280x1+250x2+400x3<60000 28 end gin 3 说明: "gin 3"表示"前3个变量为整数", 等价于: gin x1 gin x2 gin x3 求解结果 OBJECTIVE FUNCTION VALUE 1) 632.0000 VARIABLE VALUE REDUCED COST X1 64.000000 -2.000000 X2 168.000000 -3.000000 X3 0.000000 -4.000000 最优解 x1=64, x2=168, x3=0, 最优值 z=632. 若生产某一类型汽车, 则至少生产 80 辆, 求生产计划. 方法 1: 修改模型为 max z=2x1+3x2+4x3 s.t. 1.5x1+3x2+5x3≤600 280x1+250x2+400x3≤60000 x1,x2,x3 =0 或≥80 将最后一个条件分解为 8 个LP 子模型, 其中 3 个子模型应去掉, 然后逐一求解, 比较目标函数 值, 再加上整数约束, 得最优解: x1=80, x2= 150, x3=0, 最优值 z=610. x1=0, x2=0, x3≥80 x1=0, x2≥80, x3=0 x1=0, x2≥80, x3≥80 (*) x1≥80, x2=0, x3=0 x1≥80, x2≥80, x3=0 x1≥80, x2=0, x3≥80 x1≥80, x2≥80, x3≥80 (*) x1=0, x2=0, x3=0 (*) 方法 2: 引入 0-1 变量, 化为整数规划. x1=0 或≥80 变为 x1≤M·y1, x1≥80y1, y1∈{0,1} x2=0 或≥80 变为 x2≤M·y2, x2≥80y2, y2∈{0,1} x3=0 或≥80 变为 x3≤M·y3, x3≥80y3, y3∈{0,1} M 为大的正数, 可取 1000. 在LINDO 中对 0-1 变量的限定: int y1 int y2 int y3 解的结果 OBJECTIVE FUNCTION VALUE 1) 610.0000 29 VARIABLE VALUE REDUCED COST X1 80.000000 -2.000000 X2 150.000000 -3.000000 X3 0.000000 -4.000000 Y1 1.000000 0.000000 Y2 1.000000 0.000000 Y3 0.000000 0.000000 最优解同前 例2原油采购与加工 某公司用两种原油(A 和B)混合加工成两种汽油(甲和乙). 甲、乙两种汽油含原油 A 的最低比例 分别为 50%和60%, 每吨售价分别为 4800 元和 5600 元. 该公司现有原油 A 和B的库存量分别为 500 吨和 1000 吨, 还可以从市场上买到不超过 1500 吨的原油 A. 原油 A 的市场价为: 购买量不超过 500 吨时的单价为 10000 元/吨;购买量超过 500 吨但不超过 1000 吨时, 超过 500 吨的部分 8000 元/ 吨; 购买量超过 1000 吨时, 超过 1000 吨的部分 6000 元/吨. 该公司应如何安排原油的采购和加工? 问题分析 市场上可买到不超过 1500 吨的原油 A: ? 购买量不超过 500 吨时的单价为 10000 元/吨; ? 购买量超过 500 吨但不超过 1000 吨时, 超过 500 吨的部分 8000 元/吨; ? 购买量超过 1000 吨时, 超过 1000 吨的部分 6000 元/吨应如何安排原油的采购和加工? ? 利润: 销售汽油的收入- 购买原油 A 的支出. ? 难点: 原油 A 的购价与购买量的关系较复杂. 模型建立 决策变量: 设原油 A 的购买量为 x 吨, 用原油 A 生产汽油甲, 乙的数量分别为 x11, x12, 用原油 B 生产汽油甲,乙的数量分别为 x21, x22 , c(x) 是购买原油 A 的支出(千元) 目标函数: 利润(千元) max z =4.8(x11+ x21)+5.6(x12+ x22)- c(x) c(x)如何表述?由已知条件得 10 , 0 500 ( ) 8 1000, 500 1000 6 3000, 1000 1500 x x c x x x x x ≤ ≤ ? ? = + ≤ ≤ ? ? + ≤ ≤ ? 约束条件: 库存 500 吨A, 购买 x 吨A, 库存 1000 吨B. 11 12 21 22 500 1000 1500 x x x x x x + ≤ + + ≤ ≤ 30 11 11 21 11 21 12 12 22 12 22 0.5 0.6 2 3 x x x x x x x x x x ≥ ? ≥ + ≥ ? ≥ + 目标函数中 c(x)不是线性函数, 是非线性规划; 对于用分段函数定义的 c(x), 一般的非线性规划软件也难以输入和求解; 想办法将模型化简, 用现成的软件求解. 模型求解: 方法 1 x1 , x2 , x3 是以价格 10, 8, 6(千元/吨)采购 A 的吨数 x= x1+x2+x3, c(x) = 10x1+8x2+6x3 目标函数: max z =4.8(x11+ x21)+5.6(x12+ x22)- (10x1+8x2+6x3 ) 500 吨≤x≤1000 吨, 超过 500 吨的 8 千元/吨 增加约束 只有当以 10 千元/吨的价格购买 x1=500(吨)时, 才能以 8 千元/吨的价格购买 x2, 1 2 2 3 1 2 3 ( 500) 0, ( 500) 0, 0 , , 500 x x x x x x x 非线性规划模型, 可以用 LINGO 求解. Model: max= 4.8*x11 + 4.8*x21 + 5.6*x12 +5.6*x22 - 10*x1 - 8*x2 - 6*x3; x11+x12 < x + 500; x21+x22 < 1000; 0.5*x11 - 0.5*x21 > 0; 0.4*x12 - 0.6*x22 > 0; x=x1+x2+x3; (x1 - 500) * x2=0; (x2 - 500) * x3=0; x1 < 500; x2 < 500; x3 < 500; x > 0; x11 > 0; x12 > 0; x21 > 0; x22 > 0; x1 > 0; x2 > 0; x3 > 0; end 求解结果 Objective value: 4800.000 Variable Value Reduced Cost X11 500.0000 0.0000000E+00 31 X21 500.0000 0.0000000E+00 X12 0.0000000E+00 0.0000000E+00 X22 0.0000000E+00 0.0000000E+00 X1 0.1021405E-13 10.00000 X2 0.0000000E+00 8.000000 X3 0.0000000E+00 6.000000 X 0.0000000E+00 0.0000000E+00 用库存的500吨原油A、 500吨原油B生产汽油甲, 不购买新的原油A, 利润为4,800千元. LINGO 得到的是局部最优解, 还能得到更好的解吗? 方法 2 y1, y2 , y3=1——以价格 10, 8, 6(千元/吨)采购 A, x1 , x2, x3——以价格 10, 8, 6(千元/吨)采购 A 的吨数 增加约束 500y2≤x1≤500y1, 500y3≤x2≤500y2 x3≤500 y3 y1,y2,y3 是0-1 变量, y=0→x=0, x>0→y=1 这是 0-1 线性规划模型, 可用 LINDO 求解, 购买 1000 吨原油 A, 与库存的 500 吨原油 A 和1000 吨 原油 B 一起, 生产汽油乙, 利润为 5000 千元. OBJECTIVE FUNCTION VALUE 1) 5000.000 VARIABLE VALUE REDUCED COST Y1 1.000000 0.000000 Y2 1.000000 2200.000000 Y3 1.000000 1200.000000 X11 0.000000 0.800000 X21 0.000000 0.800000 X12 1500.000000 0.000000 X22 1000.000000 0.000000 X1 500.000000 0.000000 X2 500.000000 0.000000 X3 0.000000 0.400000 X 1000.000000 0.000000 优于方法 1 的结果 第四节 接力队选拔和选课策略 本节介绍一类指派问题, 指派问题一般指若干项任务分给一些候选人来完成, 每人的专长不同, 完成每项任务取得的效益或需要的资源就不同, 如何分派任务使获得的总效益最大, 或付出的总资 源最少. 若干种策略供选择, 不同的策略得到的收益或付出的成本不同, 各个策略之间有相互制约关系, 如何在满足一定条件下作出抉择, 使得收益最大或成本最小. 例1混合泳接力队的选拔 32 5 名候选人的百米成绩如下表: 甲乙丙丁戊蝶泳 1'06"8 57"2 1'18" 1'10" 1'07"4 仰泳 1'15"6 1'06" 1'07"8 1'14"2 1'11" 蛙泳 1'27" 1'06"4 1'24"6 1'09"6 1'23"8 自由泳 58"6 53" 59"4 57"2 1'02"4 如何选拔队员组成 4*100 米混合泳接力队? 丁的蛙泳成绩退步到 1'15"2; 戊的自由泳成绩进步到 57"5, 组成接力队的方案是否应该调整? 问题分析 穷举法: 组成接力队的方案共有 5!=120 种. 0-1 规划模型, cij(秒)——队员 i 第j种泳姿的百米成绩. cij i=1 i=2 i=3 i=4 i=5 j=1 66.8 57.2 78 70 67.4 j=2 75.6 66 67.8 74.2 71 j=3 87 66.4 84.6 69.6 83.8 j=4 58.6 53 59.4 57.2 62.4 若选择队员 i 参加泳姿 j 的比赛, 记xij=1, 否则记 xij=0. 4 5 1 1 min ij ij j i z c x = = = ∑∑ 每人最多入选泳姿之一, 每种泳姿有且只有 1 人45111, 1, ,5 1, 1, ,4 ij ij j i x i x j = = ≤ = = = ∑ ∑ " " 模型求解 输入 LINDO 求解(中间省略了部分行, 实际求解时是不可省略的). MIN 66.8x11+75.6x12+87x13+58.6x14 +…… +67.4x51+71 x52+83.8x53+62.4x54 SUBJECT TO x11+x12+x13+x14 <=1 …… x41+x42+x43+x44 <=1 x11+x21+x31+x41+x51 =1 …… x14+x24+x34+x44+x54 =1 END INT 20 最优解: x14 = x21 = x32 =x43 = 1, 其它变量为 0, 成绩为 253.2(秒)=4'13"2. 方案为: 甲~ 自由泳、 乙~ 蝶泳、丙~ 仰泳、丁~ 蛙泳. 讨论:丁蛙泳 c43 =69.6→75.2, 戊自由泳 c54=62.4→57.5, 方案是否调整?(敏感性分析) IP 规划一般没有与 LP 规划相类似的理论, LINDO 输出的敏感性分析结果通常是没有意义的. 将c43, c54 的新数据重新输入模型, 用LINDO 求解. 最优解: x21=x32= x43 = x51= 1, 成绩为 4'17"7. 方案为: 乙~ 蝶泳、丙~ 仰泳、丁~ 蛙泳、戊~ 自由泳. 33 例2选课策略 某专业课程设置如下表所示: 课号 课名 学分 所属类别 先修课要求 1 微积分 5 数学 2 线性代数 4 数学 3 最优化方法 4 数学; 运筹学 微积分; 线性代数 4 数据结构 3 数学; 计算机 计算机编程 5 应用统计 4 数学; 运筹学 微积分; 线性代数 6 计算机模拟 3 计算机; 运筹学 计算机编程 7 计算机编程 2 计算机 8 预测理论 2 运筹学 应用统计 9 数学实验 3 运筹学; 计算机 微积分; 线性代数 要求至少选两门数学课、 三门运筹学课和两门计算机课, 为了选修课程门数最少, 应学习哪些课 程?选修课程最少, 且学分尽量多, 应学习哪些课程? 0-1 规划模型 决策变量 xi=1 是选修课号 i 的课程( xi=0 表示不选). 目标函数 选修课程总数最少 9 1 min i i z x = = ∑ 约束条件 最少 2 门数学课 x1+x2+x3 +x4+x5≥2 3 门运筹学课 x3+x5+x6 +x8+x9≥3 2 门计算机课 x4+x6 +x7+x9≥2 先修课要求 x3=1 必有 x1= x2=1 → x3 ≤x1, x3 ≤x2 → 2x3-x1-x2 ≤0. x4 ≤x7 → x4-x7≤0. 2x5-x1-x2 ≤0 x6-x7≤0 x8-x5≤0 2x9-x1-x2 ≤0 模型求解(LINDO) 最优解: x1 = x2 = x3 = x6 = x7 = x9=1, 其它为 0, 6 门课程, 总学分 21. 讨论: 选修课程最少, 学分尽量多, 应学习哪些课程? 课程最少 9 1 min i i z x = = ∑ 学分最多 1 2 3 4 5 6 7 8 9 max 5 4 4 3 4 3 2 2 3 w x x x x x x x x x 两目标(多目标)规划——min {z,-w} . 多目标优化的处理方法: 化成单目标优化. 以课程最少为目标, 不管学分多少, 则最优解如上, 6 门课程, 总学分 21. 以学分最多为目标, 不管课程多少, 则最优解显然是选修所有 9 门课程. 34 多目标规划 在课程最少的前提下以学分最多为目标. 增加约束 9 1 6 i i x = = ∑ 以学分最多为目标求解, 最优解: x1 = x2 = x3 =x5 = x7 = x9 =1, 其它为 0, 总学分由 21 增至 22. 注意: 最优解不唯一! 可将 x9 =1 易为 x6 =1. LINDO 无法告诉优化问题的解是否唯一. 课号 课名 学分 *1* 微积分 5 *2* 线性代数 4 *3* 最优化方法 4 4 数据结构 3 *5* 应用统计 4 6 计算机模拟 3 *7* 计算机编程 2 8 预测理论 2 *9* 数学实验 3 对学分数和课程数加权形成一个目标, 如三七开. 1 2 min 0.3 0.7 Y z w z w λ λ = ? = ? 其中 9 1 i i z x = = ∑ , 1 2 3 4 5 6 7 8 9 5 4 4 3 4 3 2 2 3 w x x x x x x x x x 最优解: x1 = x2 = x3 = x4= x5 = x6 = x7 = x9 =1, 其它为 0; 总学分 28 课号 课名 学分 1* 微积分 5 2* 线性代数 4 3* 最优化方法 4 4* 数据结构 3 5* 应用统计 4 6* 计算机模拟 3 7* 计算机编程 2 8 预测理论 2 9* 数学实验 3 讨论与思考 1 2 1 2 1 2 min , 1, 0 , 1 Y z w λ λ λ λ λ λ 其中 9 1 i i z x = = ∑ , 1 2 3 4 5 6 7 8 9 5 4 4 3 4 3 2 2 3 w x x x x x x x x x 1 2 3 λ < 最优解与 1 2 0, 1 λ λ = = 的结果相同——学分最多. 1 3 4 λ < 最优解与 1 2 1, 0 λ λ = = 的结果相同——课程最少. 35 第五节 饮料厂的生产与检修 例1饮料厂的生产与检修计划 某饮料厂生产一种饮料用以满足市场需求. 该厂销售科根据市场预测, 已经确定了未来四周该 饮料的需求量. 计划科根据本厂实际情况给出了未来四周的生产能力和生产成本, 如表中所示. 每 周当饮料满足需求后有剩余时, 要支出存贮费, 为每周每千箱饮料 0.2 千元. 问应如何安排生产计划, 在满足每周市场需求的条件下, 使四周的总费用(生产成本与存贮费之和)最小? 如果工厂必须在未来四周的某一周中安排一次设备检修, 检修将占用当周 15 千箱的生产能力, 但会使检修以后每周的生产能力提高 5 千箱, 则检修应该安排在哪一周. 周次 需求量(千箱) 生产能力(千箱) 成本(千元/千箱) 1 15 30 5.0 2 25 40 5.1 3 35 45 5.4 4 25 20 5.5 合计 100 135 问题分析 ? 除第 4 周外每周的生产能力超过每周的需求; ? 生产成本逐周上升; ? 前几周应多生产一些. 模型假设 ? 饮料厂在第 1 周开始时没有库存; ? 从费用最小考虑, 第4周末不能有库存; ? 周末有库存时需支出一周的存贮费; ? 每周末的库存量等于下周初的库存量. 模型建立 决策变量 x1 ~x4: 第1~4 周的生产量. y1 ~ y3: 第1~3 周末库存量. 存贮费: 0.2 (千元/周?千箱). 目标函数 min z=5x1+5.1x2+5.4x3+5.5x4+0.2(y1+y2+y3) 约束条件 x1-y1=15 x2+ y1 -y2=25 x3+ y2 –y3=35 x4+ y3=25 x1≤30, x2≤40 x3≤45, x4≤20 x1, x2, x3 , x4, y1, y2, y3≥0 模型求解 LINDO 求解, 最优解: x1 ~x4: 15, 40, 25, 20; y1 ~ y3: 0, 15, 5 . 36 周次 需求 产量 库存 能力 成本 1 15 15 0 30 5.0 2 25 40 15 40 5.1 3 35 25 5 45 5.4 4 25 20 0 20 5.5 4 周生产计划的总费用为 528 (千元). 检修计划 在4周内安排一次设备检修, 占用当周 15 千箱生产能力, 能使检修后每周增产 5 千箱, 检修应 排在哪一周? 检修安排在任一周均可. 0-1 变量 wt: wt=1~ 检修安排在第 t 周(t=1,2,3,4) 约束条件 能力限制 x1≤30 x1+15w1≤30 x2≤40 x2+15w2≤40+5w1 x3≤45 x3+15w3≤45+5w2+5w1 x4≤20 x4+15w4≤20+5w3+5w2+5w1 目标函数不变, 增加约束条件: 检修 1 次w1+w2+w3+w4=1 LINDO 求解 最优解: w1=1, w2 , w3, w4=0; x1~ x4: 15, 45, 15, 25; y1~ y3: 0, 20, 0 . 总费用由 528 降为 527(千元) 检修所导致的生产能力提高的作用,需要更长的时间才能得到充分体现. 例2饮料的生产批量问题 饮料厂使用同一条生产线轮流生产多种饮料. 若某周开工生产某种饮料, 需支出生产准备费 8 千元. 某种饮料 4 周的需求量、生产能力和成本 周次 需求量(千箱) 生产能力(千箱) 成本(千元/千箱) 1 15 30 5.0 2 25 40 5.1 3 35 45 5.4 4 25 20 5.5 合计 100 135 存贮费: 每周每千箱饮料 0.2 (千元). 安排生产计划, 满足每周的需求, 使4周总费用最小. 符号假设 ct ~时段 t 生产费用(元/件); ht ~时段 t (末)库存费(元/件) st~时段 t 生产准备费(元); 37 dt ~时段 t 市场需求(件); Mt ~时段 t 生产能力(件). 决策变量 xt ~时段 t 生产量; yt~时段 t (末)库存量 wt =1 ~时段 t 开工生产 (wt =0 ~不开工). 假设初始库存为 0, 制订生产计划, 满足需求, 并使 T 个时段的总费用最小. 目标函数 1 min ( ) T t t t t t t t z s w c x h y = = + + ∑ 约束条件 1 0 1, 0 , 0, 0 0, , 0, 1,..., t t t t t t t t t T t t y x y d x w x M x y y x y t T ? + ? = > ? = ≤ ? = ? = = ≥ = 这是一个混合 0-1 规划模型, 用LINDO 求解 最优解: x1~ x4 : 15, 40, 45, 0; 总费用: 554.0(千元) 思考题 1. 钢管下料问题 原料钢管: 每根 19 米 客户需求: 4 米50 根6米20 根8米15 根 问题 1. 如何下料最节省? 节省的标准是什么? 问题 2. 客户增加需求: 5 米10 根 由于采用不同切割模式太多, 会增加生产和管理成本, 规定切割模式不能超过 3 种. 如何下料最 节省? 2. 发动机生产计划 某厂向用户提供发动机.合同规定, 第一﹑二﹑三季度末分别交货 40 台﹑60 台﹑80 台.每季 度的生产费用为 f(x)=ax+bx2 (元), 其中 x 是该季度生产的台数. 若交货后有剩余, 可用于下季度交货, 但需支付存储费, 每台每季度 c 元. 已知工厂每季度最大生产能力为 100 台, 第一季度开始时无存货, 设a=50, b=0.2, c=4, 问工厂应如何安排生产计划, 才能既满足合同又使总费用最低?讨论 a,b,c 变 化对计划的影响, 并作出合理的解释. 38 第三章 微分方程模型 建立微分方程模型要对研究对象作具体分析. 一般有以下三种方法: ①根据规律建模; ②用微 元法建模; ③用模拟近似法建模. 第一节 根据规律建模 在数学、力学、物理、化学等学科中已有许多经过实践检验的规律和定律, 如牛顿运动定律、 基尔霍夫电流及电压定律、物质的放射规律、曲线的切线性质等, 这些都涉及到某些函数的变化率. 我们就可以根据相应的规律, 列出常微分方程. 下面以目标跟踪问题为例介绍. 设位于坐标原点的甲舰向位于 x 轴上点 A(1,0)处的乙舰发射导弹, 导弹始终对准乙舰. 如果乙舰 以最大的速度 V0 沿平行于 y 轴的直线行驶, 导弹的速度是 5 V0, 求导弹运行的曲线. 又乙舰行驶多远 时, 导弹将它击中? 解: 设导弹轨迹为 y=y(x), 经过时间 t, 导弹位于点 P(x,y), 乙舰位于点 Q(1,V0t). 由于导弹头始终 对准乙舰, 故此时 PQ 就是曲线 y(x)在点 P 处的切线, 因此, 由01Vtyyx?′=?得0(1 ) V t x y y ′ = ? + . 又因为弧 OP 的长度为 5|AQ|, 即20015xydx V t ′ + = ∫ 所以 2 0 1 (1 ) 1 5 x x y y y dx ′ ′ ? + = + ∫ 整理得 2 1 (1 ) 1 5 x y y ′′ ′ + = + , 并有 (0) 0, (0) 0 y y′ = = , 解得 6 4 5 5 5 5 5 (1 ) (1 ) 8 12 24 y x x 当x=1 时, 5 24 y = , 即当乙舰行到 5 (1, ) 24 处被击中, 0 0 5 24 y t V V = = . y x 0 A y=y(x) P(x,y) Q(1,V0t) 39 第二节 微元法建模 微元法建模实际上是寻求一些微元之间的关系式. 与第一种方法不同之处在于这里不是直接对 未知函数及其导数应用规律和定理来求关系式, 而是对某些微元来应用规律. 以容器漏水问题为例. 有高为1米的半球形容器, 水从它的底部小孔流出. 小孔横截面为1cm2 . 开始时容器内盛满了水, 求水从小孔流出过程中容器里水面的高度 h(水面与孔口中心距离)随时间 t 变化的规律. 解: 由流体力学知识知道, 水从孔口流出的流量 Q 可用下列公式计算: 0.62 2 dv Q S gh dt = = , 其中 0.62 为流量系数, S 为孔口横截面积.现S=1cm2 , 故0.62 2 dv gh dt = . 另一方面, 设在[ ] , t t t + ? 内, 水面高度由 h 降至 ( 0) h dh dh + < , 则2dv r dh π = ? 其中 r 是时刻 t 的水面半径. 因为 2 2 2 100 (100 ) 200 r h h h 所以 2 (200 ) dv h h dh π = ? ? 于是, 2 0.62 2 (200 ) ghdt h h dh π = ? ? . 由此得 3 1 2 2 ( 200 ) 0.62 2 dt h h dh g π = ? 满足: 0 100 t h = = . 解得 3 5 2 2 5 3 (7 10 10 3 ) 4.65 2 t h h g π = * ? + 此即容器内水面高度 h 与时间 t 之间的函数关系式. 第三节 模拟近似法建模 在社会科学、生物学、医学、经济学等学科的实践中, 常常要用模拟近似法来建立微分方程模 型. 这是因为, 这些学科中的一些现象的规律我们还不是很清楚, 即使有所了解也并不全面, 因此, 要用数学模型进行研究只能在不同的假设下去模拟实际的现象. 然后再把解的结果同实际情况作对 比. 以交通管理问题为例. 在交通十字路口, 都会设置红绿灯. 为了让那些正行驶在交叉路口或离交叉路口太近而无法停 100 r h+dh dv 0 40 下的车辆通过路口, 红绿灯转换中间还要亮起一段时间的黄灯. 对于一位驶近交叉路口的驾驶员来 说, 万万不可处于这样的进退两难的境地: 要安全停车则离路口太近; 要想在红灯亮之前通过路口 又觉得太远. 那么, 黄灯应亮多长时间才最为合理呢? 解: 各段时间应该满足以下关系: 黄灯状态应持续的时间=驾驶员反应时间+车通过交叉路口时间+通过刹车距离的时间. 设0v—表示法定速度, I —交叉路口宽度, L—典型车身长度. 则通过路口的时间为 0 I L v + , (车尾 通过路口).下面计算刹车距离. 设w—为汽车的重量, μ—磨擦系数, 则摩擦力=μw, 汽车在停车过程中, 行驶距离 x 与时间 t 的关系可由下面微分方程求得 w dt x d g w ? ? = 2 2 (F=ma). 满足: 0 0 0 0, t t x x v = = ′ = = , 于是刹车距离就是直到速度 v=0 时汽车驶过的距离, 由上式得 2 0 1 2 x gt v t ? = ? + . 令0x′ = , 所以刹车所用时间 0 0 v t g ? = , 刹车距离 2 0 0 ( ) 2 v x t g ? = . 由上面得黄灯状态时间为 0 0 ( ) x t I L A T v + + = + 0 0 2 v I L T g v ? + = + + 其中 T 是驾驶员反应时间, A, v0 关系(如图)(即黄灯周期与法定速度的关系). 假设 T=1s, L=4.5m, I=9m, 另外, 我们取具有代表性的μ=0.2, 当v0=45, 60, 80km/h 时, 黄灯时间 如下表所示. v0(km/h) A(s) 经验法的值(s) 45 5.27 3 60 6.06 4 80 7.28 5 经验法的结果一律比我们预测的黄灯状态短些. 这使人想起, 许多交叉路口红绿灯的设计可能 使车辆在绿灯转为红灯时正处于交叉路口. A v0 41 第四节 微分方程模型实例 例1. 最优捕鱼策略 为了保护人类赖以生存的自然环境, 可再生资源(如渔业、林业资源)的开发必须适度. 一种合 理、简化的策略是, 在实现可持续收获的前提下, 追求最大产量或最佳效益. 考虑对某种鱼(鲥鱼)的最优捕捞策略: 假设这种鱼分 4 个年龄组: 称1龄鱼, …, 4 龄鱼. 各年龄组每条鱼的平均重量分别为 5.07, 11.55, 17.86, 22 .99(克); 各年龄组鱼的自然死亡率均为 0.8(1/年); 这种鱼为季节性集中产卵繁殖, 平均每条 4 龄鱼的产卵量为 1.109 *105 (个), 3 龄鱼的产卵量为这个数的一半, 2 龄鱼和 1 龄鱼不产卵, 产卵和孵 化期为每年的最后 4 个月; 卵孵化并成活为 1 龄鱼, 成活率(1 龄鱼条数与产卵总量 n 之比)为1.22 *1011 /(1.22 *1011 +n). 渔业管理部门规定, 每年只允许在产卵孵化期前的 8 个月内进行捕捞作业. 如果每年投入的捕 捞能力(如渔船数、下网次数等)固定不变, 这时单位时间捕捞量将与各年龄组鱼群条数成正比, 比例 系数不妨称捕捞强度系数. 通常使用 13mm 网眼的拉网, 这种网只能捕捞 3 龄鱼和 4 龄鱼, 其两个捕 捞强度系数之比为 0.42:1. 渔业上称这种方式为固定努力量捕捞. (1)建立数学模型分析如何实现可持续捕获(即每年开始捕捞时渔场中各年龄组鱼群条数不变), 并且在此前提下得到最高的年收获量(捕捞总重量). (2)某渔业公司承包这种鱼的捕捞业务 5 年, 合同要求 5 年后鱼群的生产能力不能受到太大破坏. 已知承包时各年龄组.鱼群的数量分别为: 122, 29.7, 10.1, 3.29(*109 条), 如果仍用固定努力量的捕捞 方式. 该公司应采取怎样的策略才能使总收获量最高. 问题分析 要求研究的问题是: 对某种鱼的最优捕捞策略. 1、鱼的情况 具体数据如下表: i mi(g) r(1/年) ui(个/条) 1 2 3 4 5.07 11.55 17.86 22.99 0.8 0.8 0.8 0.8 0 0 其中, i 表示 i 龄鱼, mi 表示 i 龄鱼的质量, r 表示 i 龄鱼的自然死亡率, ui 表示平均每条 i 龄鱼的产卵量. 如果每年投入的捕捞能力(如渔船数、下网次数)不变, 这时单位时间捕捞量将与 i 成正比, 比例 系数(称捕捞强度系数)为ki 使用 13mm 网眼的拉网, 这种网只能捕捞 3、4 龄鱼, 其两个捕捞强度系 数之比为 k3:k4=0.42:1, k1=k2=0. 渔业上称这种方式为固定努力量捕捞. 2、基本假设 假设Ⅰ: 一年中, 鱼的产卵是集中在 8 月底一次性完成, 捕捞工作只在前 8 个月进行. 假设Ⅱ: 各龄鱼(不包括 4 龄鱼)只在年末瞬时才长大一岁. 鱼卵在年终才孵化完毕, 成为 1 龄鱼. 这样在计算产卵量时, 3, 4 龄鱼的条数为 t=8/12. 3、应解决的问题 1)建立数学模型分析如何实现可持续捕捞(即每年开始捕捞时渔场中各年龄组鱼的条数不变), 且在此前提下得到最高的年收获量(捕捞总量). 2)某渔业公司承包这种鱼的捕捞业务 5 年, 合同要求 5 年后鱼群的生产能力不能受到太大破 坏.已知承包时各年龄组鱼群的数量分别为 x1, x2, x3, x4, 如果仍用固定努力量的捕捞方式, 该公司应 采取怎样的策略才能使总收获量最高. 42 记号和约定 xi(t): i 龄鱼在 t 时刻的数量(t 以年为单位, i=1,2,3,4); pi: i 龄鱼的捕捞量(i=3,4); M: 捕捞总质量; Q: 每年的产卵中能孵化成 1 龄鱼的数量; N: 每年的产卵量 模型的建立 由于鱼的数量随时间变化, 可视 xi(t)为连续函数, 它的变化与时间 t、自然死亡率 r、单位时间捕 捞量 ki、卵的成活率有关. 模型 I 定义单位死亡率 , 0.8 i i dx rx r dt = ? = 单位时间捕捞量 3 4 1 2 , : 0.42:1, 0 i i i dp k x k k k k dt = = = = 则捕捞时满足 ( ) i i i dx r k x dt = ? + 对各龄鱼存在以下方程(令k=k4, 则k3=0.42k) [ ] 1 1 ( ) 0.8 , 0,1 , dx t x t dt = ? ∈ [ ] 2 2 ( ) 0.8 , 0,1 , dx t x t dt = ? ∈ 3 3 3 3 ( ) 8 (0.8 0.42 ) , 0, , 12 ( ) 8 0.8 , ,1 , 12 dx t k x t dt dx t x t dt ? ? ? = ? + ∈ ? ? ? ? ? ? ? ? ? ? ? ? ? 4 4 4 4 ( ) 8 (0.8 ) , 0, , 12 ( ) 8 0.8 , ,1 , 12 dx t k x t dt dx t x t dt ? ? ? = ? + ∈ ? ? ? ? ? ? ? ? ? ? ? ? ? 由此可解得 1 1 2 2 3 3 4 4 3 3 4 4 (1) exp( 0.8) (0), (1) exp( 0.8) (0), 8 2 ( ) exp( (0.8 0.42 ) ) (0), 12 3 8 2 ( ) exp( (0.8 ) ) (0), 12 3 (1) exp( (0.8 0.28 )) (0), 2 (1) exp( (0.8 ) (0), 3 x x x x x k x x k x x k x k x x = ? = ? = ? + * = ? + * = ? + = ? + 43 收获量为 8 8 12 12 3 3 4 4 0 0 0.42 p kx t dt p kx t dt = = ∫ ∫ 模型Ⅱ 要实现持续收获, 即每年初, 各年龄组鱼群数量不变, 同时须满足以下等式(M 为捕捞总质量, N 为卵量, Q 为存活量), 由此可建立模型: max 3 3 4 4 M p m p m = + s.t. 1 2 1 3 2 4 3 4 (0) , (0) (1), (0) (1), (0) (1) (1). x Q x x x x x x x = ? ? = ? ? = ? ? = + ? 其中 5 3 4 1 (1.109 1.109 ) 10 2 N x x 11 11 1.22 10 /(1.22 10 ) Q N N 利用约束关系, 将x1(0), …, x4(0)均用 k 来表示, 从而得出 M 关于 k 的一元函数关系: 3 3 4 4 M m p m p = + , 3 3 2 0.42 (0)(1 exp( (0.8 0.42 ) )/(0.8 0.42 ) 3 p kx k k 4 4 2 (0)(1 exp( (0.8 ) )/(0.8 ) 3 p kx k k 11 6 3 (0) 0.2463 10 3.7504 10 exp(0.28 ) (2.226 exp( 2 )) 3 k x k /(2.226 exp( 2 )) 3 k + ? * , 11 4 (0) exp( 0.28 ) 10 /(9.0272 4.0598 exp( 2 )) 3 k x k 6 1.22 10 /(0.7241 0.3252 exp( 2 )) 3 k 用计算机搜索, 得出近似解, 结果为: 年最大收获量: * 11 (0) 3.88685 10 ( ) M g = * 最佳捕捞强度系数: * 17.7600( ) K g = 此时 * 11 * 10 1 2 * 10 * 7 3 4 * 11 (0) 1.193 10 , (0) 5.360 10 , (0) 2.409 10 , (0) 7.501 10 , (0) 1.193 10 . x x x x Q = * = * = * = * = * 四、模型结果及实用性讨论 当k∈[0, 31.38]时, 鱼场可实现持续捕捞, 即满足了持续捕捞条件 M≥0. 在此前提下取得了最 优的捕捞系数: 3 龄鱼的捕捞系数为 7.4592, 4 龄鱼的捕捞系数为 17.76, 年最大捕捞量为:3.88685* 1011 (克). 44 例2. 存贮模型 某厂生产若干种产品. 轮换生产时因更换设备要付准备费; 产量大于需求时因积压资金要付贮 存费. 今已知某产品的日需求量为 100 件, 生产准备费 5000 元, 贮存费每日每件 1 元. 试安排该产品 的生产计划, 即多少天生产一次(生产周期), 每次产量多少, 使总费用最小. 设该厂生产能力非常大, 即所需数量可在很短时间内产出. 要求: 不只是回答问题, 而且要建立生产周期、产量与需求量、准备费、贮存费之间的关系. 问题分析与思考 ? 日需求 100 件, 生产准备费 5000 元, 贮存费每日每件 1 元. ? 每天生产一次, 每次 100 件, 无贮存费, 准备费 5000 元, 故每天费用为 5000 元. ? 10 天生产一次, 每次 1000 件, 贮存费 900+800+…+100 =4500 元, 准备费 5000 元, 总计 9500 元. 平均每天费用为 950 元. ? 50 天生产一次, 每次 5000 件, 贮存费 4900+4800+…+100 =122500 元, 准备费 5000 元, 总计 127500 元. 平均每天费用为 2550 元. ? 周期短, 产量小→ 贮存费少, 准备费多. ? 周期长, 产量大→ 准备费少, 贮存费多. 是否存在最佳的周期和产量, 使总费用(二者之和)最小? 这是一个优化问题, 关键在建立目标函数. 显然不能用一个周期的总费用作为目标函数. 目标函数——每天总费用的平均值. 模型假设 1. 产品每天的需求量为常数 r; 2. 每次生产准备费为 c1, 每天每件产品贮存费为 c2; 3. T 天生产一次(周期为 T), 每次生产 Q 件, 且当贮存 量降到零时, Q 件产品立即生产出来(生 产时间不计); 4. 为方便起见, 时间和产量都作为连续量处理. 建模目的 设r, c1, c2 已知, 求T, Q, 使每天总费用的平均值最小. 模型建立 离散问题连续化 将贮存量表示为时间的函数 q(t), t=0 生产 Q 件, 贮存量 q(0)=Q, q(t) 以需求 r 的速率递减, 直到 q(T)=0. 由此得 Q=rT 一周期贮存费 A c dt t q c T 2 0 2 ) ( = ∫ 一周期总费用 45 2 2 2 2 1 2 1 rT c c T Q c c C + = + = 每天总费用平均值(目标函数) 2 ) ( 2 1 rT c T c T C T C + = = 模型求解 求T使min 2 ) ( 2 1 → + = rT c T c T C . 由0dC dT = 有212rc c T = , 2 1 2 c r c rT Q = = . 模型分析 1 , c T Q ↑? ↑ , 2 , c T Q ↑? ↓ , , r T Q ↑? ↓ ↑ . 模型应用 经济批量订货公式(EOQ 公式) 用于订货、供应、存贮情形 每天需求量 r, 每次订货费为 c1,每天每件贮存费为 c2 , T 天订货一次(周期 T), 每次订货 Q 件, 且当贮存量降到 零时, Q 件立即到货. 2 1 2 rc c T = , 2 1 2 c r c rT Q = = 这是不允许缺货的存贮模型 问: 为什么不考虑生产费用?在什么条件下才不考虑? 允许缺货的存贮模型 当贮存量降到零时仍有需求 r, 出现缺货, 造成损失. 原模型假设: 贮存量降到零时 Q 件立即生产出来(或立即到货). 现假设: 允许缺货, 每天每件缺货损失费 c3 , 缺货需补足. 周期 T, t=T1 贮存量降到零. T 一周期贮存费 A c dt t q c T 2 0 2 1 ) ( = ∫ 一周期缺货费 B c dt t q c T T 3 3 1 ) ( = ∫ 46 一周期总费用 2 ) ( 2 2 1 3 1 2 1 T T r c QT c c C ? + + = 每天总费用平均值(目标函数) rT Q rT c rT Q c T c T C Q T C 2 ) ( 2 ) , ( 2 3 2 2 1 ? + + = = 求T,Q 使(,)min C T Q → 0 , 0 = ? ? = ? ? Q C T C ,为与不允许缺货的存贮模型相比, T 记作T′ , Q 记作Q′ . 3 3 2 2 1 2 ' c c c rc c T + = , 3 2 3 2 1 2 ' c c c c r c Q + = 例3. 传染病模型 随着卫生设施的改善、医疗水平的提高以及人类文明的不断发展, 诸如霍乱、天花等曾经肆虐 全球的传染性疾病已经得到有效的控制. 但是一些新的、不断变异着的传染病毒却悄悄向人类袭来. 20 世纪 80 年代十分险恶的爱滋病毒开始肆虐全球, 至今仍在蔓延;2003 年春来历不明的 SARS 病 毒突袭人间, 给人们的生命财产带来极大的危害. 长期以来, 建立传染病的数学模型来描述传染病 的传播过程, 分析受感染人数的变化规律, 探索制止传染病蔓延的手段等, 一直是各国有关专家和 官员关注的课题. 试建立传染病的数学模型, 以分析传染病的传播规律. 模型 1 已感染人数(病人) 假设每个病人每天有效接触(足以使人致病)人数为 λ, 建立如下数学模型 t t i t i t t i ? = ? ? + ) ( ) ( ) ( λ 由0)(,itiidt di = = λ 推出 t e i t i λ 0 ) ( = , ∞ → ? ∞ → i t . 若有效接触的是病人, 则不能使病人数增加, 必须区分已感染者(病人)和未感染者(健康人). 模型 2 区分已感染者(病人)和未感染者(健康人) (SI 模型) 假设: 1)总人数 N 不变, 病人和健康人的比例分别为 i(t), s(t). 2)每个病人每天有效接触人数为 λ, 且使接触的健康人致病, λ 为日接触率 建立数学模型 t t i t Ns t i t t i N ? = ? ? + ) ( )] ( [ )] ( ) ( [ λ 由si dt di λ = , 1 ) ( ) ( = + t i t s 得?????=?=0)0()1(iiiidt di λ 称为 Logistic 模型. 可解出 t e i t i λ ? ? ? ? ? ? ? ? ? ? + = 1 1 1 1 ) ( 0 47 令?????????=?11ln 0 1 i tm λ , tm 是传染病高潮到来时刻. λ (日接触率)↓ → tm↑. 1 → ? ∞ → i t , 说明病人可以治愈! 模型 3 传染病无免疫性——病人治愈成为健康人, 健康人可再次被感染(SIS 模型) 增加假设 3)病人每天治愈的比例为 ? , 即日治愈率. 建立数学模型 N i t t i t Ns t i t t Ni t t λ ? 0 (1 ) (0) di i i i dt i i λ ? ? = ? ? ? ? ? = ? λ 是日接触率, 1/? 是感染期 σ=λ/? 是一个感染期内每个病人的有效接触人数, 称为接触数. )] 1 1 ( [ σ λ ? ? ? = i i dt di , ? ? ? ? ? ≤ > ? = ∞ 1 , 0 1 , 1 1 ) ( σ σ σ i 感染期内有效接触感染的健康者人数不超过病人数. 思考 模型 2(SI 模型)如何看作模型 3(SIS 模型)的特例 模型 4 传染病有免疫性——病人治愈后即移出感染系统, 称移出者 SIR 模型 假设 1)总人数 N 不变, 病人、健康人和移出者的比例分别为 i(t), s(t), r(t) 2)病人的日接触率 λ , 日治愈率 ?, 接触数 s=λ/?. 建立数学模型 1 ) ( ) ( ) ( = + + t r t i t s 需建立 i(t), s(t), r(t)的两个方程 t t Ni t t i t Ns t i t t i N ? ? ? = ? ? + ) ( ) ( ) ( )] ( ) ( [ ? λ t t i t Ns t s t t s N ? ? = ? ? + ) ( ) ( )] ( ) ( [ λ 48 ? ? ? ? ? ? ? ? ? = = ? = ? = 0 0 ) 0 ( , ) 0 ( s s i i si dt ds i si dt di λ ? λ i0+s0≈1(通常 r(0)=r0 很小) 无法求出 i(t), s(t)的解析解, 在相平面 s-i 上研究解的性质. 由上模型得 ? ? ? ? ? = ? = = 0 0 1 1 i i s ds di s s σ 相轨线 0 0 0 ln 1 ) ( ) ( s s s i s s i σ + ? + = 定义域 } 1 , 0 , 0 ) , {( ≤ + ≥ ≥ = i s i s i s D , 在D内作相轨线 i(s)的图形, 进行分析 0, ( 1/ ) , m i s s t i i s i s σ ∞ ∞ ↓ = = = 图形 满足 0 0 0 1 ln 0 s s i s s σ ∞ ∞ 0 1 1/ s P i t σ > → 先升后降至 0, 传染病蔓延. 0 2 1/ s P i t σ < → 单调降至 0, 传染病不蔓延. 预防传染病蔓延的手段 传染病不蔓延的条件—— 0 1/ s σ < . λ (日接触率)↓ ? 卫生水平↑. ?(日治愈率)↑ ? 医疗水平↑. 降低 0 0 0 0 0 ( 1) s s i r r 即群体免疫. σ 的估计: 由0001ln 0 s s i s s σ ∞ ∞ 忽略 i0, 得00ln ln s s s s σ ∞ ∞ ? = ? . 49 第四章 图论模型 图论最早起源于一些数学游戏的难题研究, 如欧拉所解决的哥尼斯堡七桥问题, 以及在民间广 泛流传的一些游戏难题, 如迷宫问题、 博奕问题、 棋盘上马的行走路线问题等. 这些古老的难题, 当 时吸引了很多学者的注意.在这些问题研究的基础上又继续提出了著名的四色猜想和汉米尔顿(环游 世界)数学难题. 1847 年, 图论应用于分析电路网络, 这是它最早应用于工程科学, 以后随着科学的发展, 图论 在解决运筹学, 网络理论, 信息论, 控制论, 博奕论以及计算机科学等各个领域的问题时, 发挥出越 来越大的作用.在人们的社会实践中, 图论已成为解决自然科学、工程技术、社会科学、生物技术 以及经济、 军事等领域中许多问题的有力工具之一, 因此越来越受到数学家和实际工作者的喜爱. 我 们所学的这一章只是介绍一些基本概念、原理以及一些典型的应用实例, 目的是在今后对工程技术 有关学科的学习研究时, 可以把图论的基本知识、方法作为工具. 第一节 图论的基本概念 一、图的定义 图G是指由非空有限集合 V(G) 和V(G)中某些元素的无序对的集合 E(G)构成的二员组(V(G), E(G)), 其中 V(G) 成为 G 的顶点集合, 其中的元素成为 G 的顶点. E(G) 成为 G 的边集, 其中的元 素称为 G 的边. 一般简单记为 G=(V(G), E(G)), 在不发生混淆的情况下, 可记为 G=(V, E). 如果 V={v1, …, vn}, 那么 E 中的元素 e 与V中某两个元素构成的无序对{vi, vj}相对应, 记为 e=vivj 或e=vjvi. 图可用图形来表表示, 如: 简单图: 既没有环也没有重边的图. 开始我们都是讲简单图. 二、点与点, 点与边, 边与边间的关系 若e=vivj∈E, 则称顶点是 vi 和vj 相邻, 并称 vi, vj 为边 e 的端点, 也称 e 与vi, vj 关联. 若e1, e2∈E, 且e1 和e2 有公共的端点, 则称 e1 与e2 是相邻的. 三、阶与度 1、图的阶: 一个图的顶点数. 2、 顶点的度 dG(v): 图G中和 V 关联的边的数目(与V关联的每个环算作两条边). 顶点有偶点和 奇点之分. 当dG(v)=0 时称其为孤立点. 结论 1: 设G=(V, E)是一个图, 则有 | | 2 ) ( E v d V v G = ∑ ∈ . 即奇点的个数必为偶数. 四、点边序列 若序列: v0e1v1…ekvk, 其中 ei∈E, vj∈V, i=1,2,…,k, j=0,1,…,k, 且相邻两项有关联关系, 则称此序 列是一条道路(途径、通道), 其中称 v0 为起点, 称vk 为终点, 而k称为此道路的长. 1、迹: 如果 G 中的一条道路 W 上的边互不相同, 则称 W 为G的迹. 2、链: 如果 G 中的一条道路 W 上的点互不相同, 则称 W 为G的链(路). 若W的起点、终点分 50 别为 vi, vj, 则W可简记为(vi, vj)或P(vi, vj). 3、闭通道: 如果通道的长至少为 1, 且起点和终点重合, 则称该通道为闭通道. 4、闭迹: 起点和终点重合且长至少 1 的迹称为闭迹. 5、圈: 起点和内部顶点互不相同的闭迹称为圈. 五、子图 1、子图: 如果图 G 和图 H 满足 V(H) ? V(G), E(H) ? E(G), 则称 H 是G的子图, 记为 H ? G. 2、支撑子图: 如果 V(H)=V(G), 且E(H) ? E(G), 则称 H 为G的支撑子图. 3、导出子图: 设V'是图 G 的顶点集 V(G)的一个非空子集, 令},,), ( | { V v v v v e G E e e E j i j i ′ ∈ = ′ ∈ ′ ′ = ′ 则称 G 的子图 G'=(V',E')为由 V 导出的子图, 记为 G'=G[V']. 设E' ? E(G), E'≠?令{|}VvvE′′=是中某条边的端点 则称 G 的子图 ( , ) G V E ′ ′ ′ = 为由 E'导出的子图, 记为 G'=G[E']. 4、生成子图: 设G与H是两个图, 且 V H V G E H E G = ? , 则称 H 为G的生成子图. 六、图的连通性 1、 连通图: 如果对于图 G 的任意两个顶点 vi 和vj, G 中都存在链 P(vi,vj), 则称图 G 是连通图. 否 则称为非连通图. 2、连通子图: 如果 H G ? , 且H是连通图, 则称 H 是G的连通子图. 3、极大连通子图: 指H为G的连通子图, 并且 G 中不存在连通子图 H'使得 , H H H H ′ ′ ? ≠ . 此时也称它为 G 的连通分支. 从而可知, 连通图恰有一个连通分支, 而非连通图则有两个或两个以 上的连通分支. 无向图的连通性的图例 无向图 G 无向图 G 的三个连通分量 七、连通图的应用 8 个城市 v0, v1 , …, v7 之间有一个公路网(如图所示), 每条公路为图中的边, 边上的权数表示 51 通过该公路所需的时间.设你处在城市 v0, 那么从 v0 到其他各城市, 应选择什么路径使所需的时间 最短? 本例比较简单, 可以用穷举法选择出从 v0 到其他各城市之间的最短路, 即运行时间最短的 路.然而当公路网比较复杂时, 是否有好的算法呢? 狄克斯特拉在 1959 年提出一种最短路的算法, 至今仍被认为是最好的算法之一. 设(v0, u1, u2, …, un)为从 v0 到un 的最短路, 记001121 n n n d v u v u u u u u ω ω ω ? 称之为从 v0 到un 的最短路的权数.若S为V的任意子集, 则称 ) , ( min ) , ( 0 0 u v d S v d S u∈ = 为从v0 到S的最短路的权数. 对于例1, 设S为节点集 V 的一个节点子集, v0 S S V S S c 为 ,设 \ = ∈ 的节点余集.容易看出, 若( v0, u1, u2, …, u v ? ? , ) 为从 v0 到cS的最短路, 则必有 , , c S v S u ∈ ∈ ? ? 使(v0, u1, u2, · · ·, u ? ) 为从 v0 到u ? 的最短路.于是 { } 0 0 0 0 0 0 ( , ) min c c c d v v d v S d v v d v u u v d v S d v u u v u S v S ω ω ? ? ? ? ? ? = ? ? ? = + ? ? = + ? ∈ ? ∈ ? 上式为狄克斯特拉算法的理论基础.下面以例1为例说明算法步骤: 设S0={v0}, { } 0 1 2 7 , , , c S v v v = ??? , 求{}0000000(,)min min{ ( , )} c c c d v S d v u u v d v v u S v S v S ω = + = ∈ ∈ ∈ 容易看出, 0 0 0 1 1 c d v S v v ω = = . 把v1 加入 S0, 设S1={v0, v1}, 1 2 3 7 c S v v v = ??? , 再求 0 1 0 0 2 1 1 ( , ) min{2 c c d v S d v u u v d v v u S v S ω = + = = ∈ ∈ . 52 把v2 加入 S1, 设S2={v0, v1, v2}, 2 3 4 7 c S v v v = ??? , 求02002230322(,)min 3 c c d v S d v u u v d v v v v d v v u S v S ω ω ∈ ∈ 以此类推, 可求出从 v0 到其它任意节点的最短路和它们的权数. 若把从 v0 到vj(j=1, 2, …, 7)的最 短路的权数标在该点处, 并令 v0 处的权数为零, 则所得结果如图所示. 实际上, 它是原图的以 v0 为根的生成树, 直接表明从 v0 到其它各节点的最短路径. 八、图的运算 1、减: 若E1 ? E(G), 从图 G 中删去 E1 的所有边得到得图, 记为 G-E1. 若V1 ? V(G), 且V1≠V(G), 则把从图 G 中删去 V1 得所有顶点以及 V1 中顶点关联得边得到得图 记为 G-V1. 易得 G[V']=G-(V(G)\V'). 2、 加: 若H?G, E2 ? E(G)\E(H), 且E2 中每条边得端点都属于 V(G), 在图 H 中添上 E2 得所有边 得到得图记为 H+E2. 3、并: 若H1 ? G, H2 ? G, 则把顶点集为 V(H1)∪V(H2), 边集为 E(H1)∪E(H2)的图成为 H1 和H2 的并, 记为 H1∪H2. 若E(H1)∩E(H2)=?, 则称 H1∪H2 为H1 与H2 的和, 记为 H1+H2. 九、边割 设G=(V,E)是一个图, S 为V的一个非空真子集, \ S V S = , 记 i j i j S S v v E v S v S = ∈ ∈ ∈ 如果[ , ] S S ≠ ? , 则称[ , ] S S 是G的一个边割. 若E'是G的一个边割, 但E'的任何真子集都不 是G的边割, 则称 E 为G的极小边割, 此时称之为 G 的割集. 结论 1: 设C和[ , ] S S 分别时图 G 的圈和边割, 则0(mod 2) E C S S ≡ ∩ 十、有向图的基本概念 1、有向图 D: 指由非空有限集合 V(D) 和V(D)中某些元素的无序对的集合 A(D)构成的二员组 (V(D), A(D)), 其中 V(D)成为 D 的顶点集合, 其中的元素成为 D 的顶点. A(D)成为 D 的弧集, 其中的 元素称为 D 的弧. 在不发生混淆的情况下, 可记为 D(V, A). 如果 1 2 n V D v v v = , 那么 A 中的元素 a 与V(D)中某两个元素构成的有序对{ , } i j v v 相对应, 记为 a={vi,vj}. 其中 vj 成为 a 的尾, vi 称为 a 的头, a 称为 vi 的出弧, 也称为 vj 的入弧. 53 顶点 v 的在 D 中的出弧的数目记为 ( ) D d D + , 称之为 v 的出度; 顶点 v 在D中入弧的数目记为 ( ) D d D ? , 称之为 v 的入度. 结论 3: 设D=(V,A)是一个有向图, 则 D D v V v V d v d v A + ? ∈ ∈ = = ∑ ∑ 注意: 在无向图中的定义(例如通道、迹等)也可相应的扩充为有向图的情况. 十一、几类重要的图 1、完全图: 任何一对顶点都相邻的简单图. N 阶完全图记为 Kn. 划分: 设S是一个集合, 如果存在 k 个集合 S1, S2, …, Sk, 使得 1 k i i S S = = ∪ 且()ijSSij=?≠∩则称 S1, S2, …, Sk 为S的一个划分. 2、二部图(二分图): 设G是一个图, 如果存在 V(G)的一个划分 X,Y, 使得 G 的任何一条边的一 个端点在 X 中, 另一个在 Y 中, 则称 G 为二部图, 记作 G=(X,Y,E). 完全二部图是指这样的简单二部图 G=(X,Y,E), 其中 X 的每个顶点都与 Y 的每个顶点相邻. 若||,| | X m Y n = = , 则完全二部图 G=(X,Y,E)记为Kmn. 结论 4: 图G是二部图当且仅当 G 中不包含奇圈. 3、Euler 图:是指图 G 中存在包含一切边的闭迹 W. 结论 5: 图G是Euler 图当且仅当 G 中无奇点. 4、Hamlton 图: 是指图 G 中存在一条包含所有顶点的圈 C. 结论 6: G 是一个(p,q)图, p≥3. 若,,uvV?∈有()()GGdudvp+≥,则称图 G 是Hamlton 图. 5、 无向网络(赋权图): 在边加权 w 的图 G.其中 w 是E(G)上的一个实值函数, 称为 G 的权函数. 这时G记为 G=(V(G),E(G),w). 第二节 图的存储结构 邻接矩阵和加权邻接矩阵(labeled adjacency matrix) ?无权值的有向图的邻接矩阵 设有向图具有 n 个结点, 则用 n 行n列的布尔矩阵 A 表示该有向图; 并且 A[i,j]=1,如果 i 至j有一条有向边; A[i,j]=0 如果 i 至j没有一条有向边注意: A[i,i]=0. 出度:i 行之和. 入度: j 列之和. 54 ?无权值的无向图的邻接矩阵 设无向图具有 n 个结点, 则用 n 行n列的布尔矩阵 A 表示该无向图; 并且 A[i,j]=1,如果 i 至j有一条无向边; A[i,j]=0 如果 i 至j没有一条无向边 注意: A[i,i]=0. i 结点的度:i 行或 i 列之和. 上三角矩阵或下三角矩阵. ?有向图的加权邻接矩阵 设有向图具有 n 个结点, 则用 n 行n列的矩阵 A 表示该有向图; 并且 A[i,j]=a,如果 i 至j有一条有向边且它的权值为 a. A[i,j]=空或其它标志; 如果 i 至j没有一条 有向边. 优点: 判断任意两点之间是否有边方便, 仅耗费 O(1)时间. 缺点: 即使<0), 因此每条边赋予向前费用权 ω+ (i,j)与向后费用权 ω- (i,j): ( , ) a i j f i j c i j i j f i j c i j ω+ < ? = ? +∞ = ? 若若61 0 ( , ) , ( , ) 0 a i j f i j i j f i j ω? ? > ? = ? +∞ = ? 若若对赋权后的有向图, 如果把权 ω(i,j)看作长度, 即可确定从 s 到s′的费用最小的非饱和路, 它等价于 确定从 s 到s′的最短路.确定了非饱和路后就可确定该路的最大可增流量, 因此需对每一条边确定 一个向前可增流量 ?+ ( , ) i j 与向后可增流量 ?? ( , ) i j : ( , ) 0, c i j f i j f i j c i j i j f i j c i j + ? < ? ? = ? = ? 若若0(,)0, ( , ) 0 f i j f i j i j f i j ? > ? ? = ? = ? 若若因此, 确定最小费用最大流的算法如下: (1)从零流开始, 令f0≡0; 赋权: 当),(),(jicjifk < , ? ? ? ? ? ? = ? = + + ) , ( ) , ( ) , ( ) , ( ) , ( j i f j i c j i j i a j i k ω ; 当?????=?+∞ = = + + 0 ) , ( ) , ( ), , ( ) , ( j i j i j i c j i fk ω ; 当;),(),(),(),(,0),(?????=??=>??jifjijiajijifkkω当(,)(,)0, ( , ) 0 k i j f i j i j ω? ? ? = +∞ ? = ? ? = ? ? (3)确定一条从 s 到s′ 的最短路 )} , ( , ), , ( ), , {( ) , ( 2 1 1 s i i i i s s s R k ′ ? ? ? = ′ . 若),(ssR′的长度为+∞, 则已 得最小费用最大流, 停机, 否则转入(4). (4)确定沿着该路 ) , ( s s R ′ 的最大可增流量 )} , ( , ), , ( ), , ( min{ 2 1 1 s i i i i s a k ′ ? ? ? ? ? ? = , 其中根据边的 取向决定取 ?+ 或?? . (5)生成新的流 ? ? ? ? + = + 为向后边 为向前边 ) , ( , ) , ( ) , ( , ) , ( ) , ( 1 j i a j i f j i a j i f j i f k k k 若fk+1(i,j)已为最小费用最大流, 停机, 否则转入(2). 下面举例说明. 例 求如图 4-5 所示的最小费用最大流.图中, 边旁的两个数字分别为 c(i,j)和a(i,j). 62 从零开始, f(i,j)≡0.图4-5 即可作为第一次迭代的赋权图(略去取 +∞ 的权). 第一次迭代: (1)由最短路法, 可得从 s 到s′ 的最短路为 )}, 6 , 5 ( ), 5 , 3 ( ), 3 , 1 {( ) , ( = ′ s s R 其费用为 4. (2)确定沿 ) , ( s s R ′ 的最大可增流量为 3 } 3 , 3 , 4 min{ )} 6 , 5 ( ), 5 , 3 ( ), 3 , 1 ( min{ = = ? ? ? = a . (3)调整流量, 生成新的流 3 ) 6 , 5 ( , 3 ) 5 , 3 ( , 3 ) 3 , 1 ( = = = f f f . (4)赋权, 得新的赋权图, 如图 4-6 所示.其中, 边旁的两个数字分别为(i,j)和ω(i,j). 第二次迭代: (1)从s到s′ 的最短路为 ( , ) {(1,3),(3,4),(4,6)}, R s s′ = 其中费用为 5. (2)沿),(ssR′的最大可增流量为 min{ (1,3), (3,4), (4,6)} min{1,5,7} 1 a (3)调整流量, 生成新的流 f(1,3)=4, f(3,4)=1, f(4,6)=1. (4)赋权, 得新的赋权图, 如图 4-7 所示. 第三次迭代: (1)从s到s′的最短路为 R(s, s′)={(1,2),(2,5),(5,3),(3,4),(4,6)}, 费用为 6. (2) ) , ( s s R ′ 的最大可增流量 3 } 6 , 4 , 3 , 4 , 3 min{ )} 6 , 4 ( ), 4 , 3 ( ), 3 , 5 ( ), 5 , 2 ( ), 2 , 1 ( min{ = = ? ? ? ? ? = a . (3)调整流量, 生成新的流 f(1,2)=3, f(2,5)=3, f(5,3)=0, f(3,4)=4, f(4,6)=4. (4)赋权得新的赋权图, 如图 4-8 所示. 至此, 在图 4-8 中, 不再存在从 s 到s′的最短路, 因此得最小费用的最大流为 f(1,2)=3, f(2,5)=3, f(3,4)=4, f(4,6)=4, f(1,3)=4, f(5,6)=3, f(2,4)=0, f(3,5)=0. 最小费用∑f(i,j)a(i,j)=35. 图4-5 图4-6 图4-7 图4-8 2. 运输问题 某产品有 m 个产地, 第i个产地的月最大产量和单位产品的成本分别为 ci 单位/月和 ai(i=1,2,…,m).该产品有 n 个需求地, 分别的需求量为 bj 单位/月(j=1,2,…,n).又已知从第 i 个产地到 第j个需求地的最大可能运输量为 ω(i,j)单位/月, 单位运价为 l(i,j).问各地应生产多少产品?如何调 63 运, 才能尽可能满足需求, 并使其总开支最省? 首先, 我们虚设一个发点 s 和收点 s′ .从发点到第 i 个产地的权数设为 ci, 单位费用为 ai(i=1,2,…,m). 从产地 i 到需求地 j 的权数 ω(i,j)其单位费用为 l(i,j) (i=1,2,…,m , j=1,2,…,n), 从需求地 j 到收点 s′ 的权数为 bj, 单位成本为零. 这样就化为图 4-9 的网络的最小费用最大流问题. 运输问题化为最小费用最大流问题的技巧在于虚设一个发点和收点, 以及确定各边相应的权数 与单位费用. 另外, 上述运输问题也可表述为线性规划问题, 即下面的总费用最小的问题: 1 1 1 1 min 1,2, , ( , ) , 1,2, , ( , ) 0, 1,2, , , 1,2, , m n i j n i j m j i l i j g i j st g i j c i m g i j b j n g i j i m j n = = = = ? ? ? ? ≤ = ??? ? ? ? ? ? ? ∑∑ ∑ ∑ 第一、二个约束条件蕴含着 1 1 1 1 ( , ) n m n m j i j i j i b g i j c = = = = ≤ ≤ ∑ ∑∑ ∑ . 即总需求量小于等于总供应量. 第五节 图的遍历 简单可这样说, 从连通图的一顶点出发每边通过一次或每个顶点恰好通过一次. 并不是每个连 通图上都能实现的事; 称图中含每条边的行迹为 Euler 行迹, 闭的 Euler 行迹叫做 Euler 回路, 含Euler 回路的图叫做 Euler 图. 图的遍历有两种基本方法: 1、 深度优先遍历 从图中指定的起点 v0 出发(先访问 v0)访问它的任意相邻接的顶点 v1, 再访问 v1 的任一个未访问 的相邻接顶点 v2, 如此下去, 直到某顶点已无被访问过的顶点, 则返回一步找前一个顶点的其他沿未 被访问的邻接顶点, 重复以上过程直到所有的顶点都被访问. 例如, 按深度优先遍历法遍历下图, 遍图4-9 64 历的结果为: v0 v1 v2 v5 v4 v6 v3 v7. 2、 广度优先遍历 从图指定的起点 v0 出发, 访问与它相邻的所有顶点 v1, v2,… , 然后再依次访问 v1, v2,… , 邻接的 尚未被访问的所有顶点, 再从这些顶点出发访问与它们相邻接的尚未被访问的顶点, 直到所有顶点 均被访问过为止. 例如, 按广度优先遍历法遍历下图, 遍历的结果为: v0 v1 v3 v4 v2 v6 v5 v7. 关于 Euler 遍历问题有下列的结论: (1)关于图是 Euler 图的充要条件. (2)连通图 G 有Euler 行迹的充要条件是连通图 G 中至少有两个奇次点. Euler 图解决了"七桥问题": 从河的一岸或一个岛上出法, 不重复地走遍七座桥. 这个问题的一个拓展问题就是"中国邮递员"问题: 邮递员每天要走过固定的区域内的所有街 道. 所谓欧拉回路与哥尼斯堡 7 桥问题相联系的.在哥尼斯堡 7 桥问题中, 欧拉证明了不存在这样 的回路, 使它经过图中每条边且只经过一次又回到起始点. 与此相反, 设G(V, E)为一个图, 若存在一 条回路, 使它经过图中每条边且只经过一次又回到起始点, 就称这种回路为欧拉回路, 并称图 G 为 欧拉图. 在一个图中, 连接一个节点的边数称为该节点的度数.对欧拉图, 我们有下列结果: 定理 1 对连通图 G(V, E), 下列条件是相互等价的: 65 (1)G 是一个欧拉图; (2)G 的每一个节点的度数都是偶数; (3)G 的边集合 E 可以分解为若干个回路的并. 证明 : (1)→(2)已知 G 为欧拉图, 则必存在一个欧拉回路.回路中的节点都是偶度数. (2)→(3)设G中每一个节点的度数均为偶数. 若能找到一个回路 C1 使G=C1, 则结论成立. 否则, 令G1=G\C1, 由C1 上每个节点的度数均为偶数, 则G1 中的每个节点的度数亦均为偶数. 于是在 G1 必存 在另一个回路 C2.令G2=G1\C2, …., 由于 G 为有限图, 上述过程经过有限步, 最后必得一个回路 Cr 使Gr=Gr-1\Cr 上各节点的度数均为零, 即Cr=Gr-1.这样就得到 G 的一个分解 G=C1∪C2∪…∪Cr. (3)→(1)设G=C1∪C2∪…∪Cr, 其中 Ci=(i=1,2,…,r)均为回路.由于 G 为连通图, 对任意回路 Ci, 必存在另一个回路 Cj 与之相连, 即Ci 与Cj 存在共同的节点. 现在我们从图 G 的任意节点出发, 沿着 所在的回路走, 每走到一个共同的节点处, 就转向另一个回路, …, 这样一直走下去就可走遍 G 的每 条边且只走过一次,最后回到原出发节点, 即G为一个欧拉图. 利用定理 1 去判断一个连通图是否为欧拉图比较容易, 但要找出欧拉回路, 当连通图比较复杂 时就不太容易了.下面介绍一种求欧拉回路的算法. 弗罗莱算法 算法步骤如下: (1)任取起始点 v0, v0→R. (2)设路 1 1 2 1 1 0 2 r r i i i r i i R e v v e v v e v v ? = ??? 已选出, 则从 E\{e1, e2, …er,}中选出边 er+1, 使er+1 与vi 相连, 除非没有其它选择, Gr\{ er+1}仍应为连通的. (3)重复步骤(2), 直至不能进行下去为止. 定理 2 连通的有向图存在欧拉回路的充分必要条件是对任意节点, 进入该节点边数(进数)与 离开该点的边数(出数)相等. 66 第五章 数学软件使用初步 在历届的数学建模竞赛中, 参赛者对数学软件的熟练掌握程度对建模的成功与失败起着非常重 要的作用. 从数据的处理到公式的推导, 从模型的建立到模型的求解, 没有数学软件的帮助几乎是 不可能完成的. 从今后的趋势来看, 数学软件在建模过程中会占据更为重要的地位. 现在数学软件 能做的事情已经非常之多, 几乎涉及到了所有的数学领域, 数的运算、 公式推导、 方程求解、 微积分、 线性代数、概率统计、规划求解、函数作图、图论、数论、数值计算、信号处理等等. 随着计算技 术的不断发展, 数学软件的功能会更强大、更智能. 当今, 综合性的数学软件主要有 MATLAB, MAPLE, Mathematica, MathCAD 等, 专用数学软件 主要有 SPSS, LINDO, LINGO, 这些在建模竞赛过程中都可能用到. 在这里, 我们主要介绍 MATLAB、MAPLE、LINDO 和LINGO. 第一节 MATLAB §1.1 MATLAB 简介 MATLAB 名字是由 MATrix 和LABoratory 两个词的前三个字母组合而成的. 它是 MathWorks 公 司于 1982 年推出的一套高性能的数值计算和可视化数学软件. 由于使用编程运算与人进行科学计算 的思路和表达方式完全一致, 所以不象学习其它高级语言如 Basic、fortran 和C等那样难于掌握. 用Matlab 编写程序犹如在演算纸上排列出公式与求解问题, 所以又被称为演算纸式科学算法语言. 在 这个环境下, 对所要求解的问题, 用户只需简单地列出数学表达式, 其结果便以数值或图形方式显 示出来. MATLAB 的含义是矩阵实验室, 主要用于方便矩阵的存取, 其基本元素是无须定义维数的矩阵. 自问世以来, 就以数值计算称雄. MATLAB 进行数值计算的基本单位是复数数组(或称阵列), 这使得 MATLAB 高度 "向量化" . 经过十几年的完善和扩充, 现已发展成为线性代数课程的标准工具. 由于 它不需定义数组的维数, 并给出矩阵函数、 特殊矩阵专门的库函数, 使之在求解诸如信号处理、 建模、 系统识别、控制、优化等领域的问题时, 显得大为简捷、高效、方便, 这是其它高级语言所不能比拟 的. 美国许多大学的实验室都安装有供学习和研究之用. 在那里, MATLAB 是攻读学位的大学生、硕 士生、博士生必须掌握的基本工具. MATLAB 中包括了被称作工具箱(TOOLBOX)的各类应用问题的求解工具. 工具箱实际上是对 MATLAB 进行扩展应用的一系列 MATLAB 函数(称为 M 文件), 它可用来求解各类学科的问题, 包 括信号处理、图象处理、控制系统辨识、神经网络等. 随着 MATLAB 版本的不断升级, 其所含的工 具箱的功能也越来越丰富, 因此, 应用范围也越来越广泛, 成为涉及数值分析的各类工程师不可不 用的工具. §1.2 基本运算与函数 在MATLAB 下进行基本数学运算, 只需将运算式直接打入提示号(>>)之后, 并按入 Enter 键即可. 例如: >>(5*2+1.3-0.8)*10/25 ans = 67 4.2000 MATLAB 会将运算结果直接存入一变量 ans, 代表 MATLAB 运算后的答案(Answer), 并显示其 数值于屏幕上. (为简便起见, 在下述各例中, 我们不再写出 MATLAB 的提示号. ) 小提示 ">>"是MATLAB 的提示符号(Prompt), 但在 PC 中文视窗系统下, 由于编码方式不同, 此提示符 号常会消失不见, 但这并不会影响到 MATLAB 的运算结果. 我们也可将上述运算式的结果设定给另一个变量 x: x = (5*2+1.3-0.8)*10^2/25 x = 42 此时 MATLAB 会直接显示 x 的值. 由上例可知, MATLAB 认识所有一般常用到的加(+)、减(-)、 乘(*)、除(/)的数学运算符号, 以及幂次运算(^). 小提示 MATLAB 在使用变量之前可以不事先经过变量声明(Variable declaration). 若不想让 MATLAB 每次都显示运算结果, 只需在运算式最后加上分号";"即可, 如下例: y = sin(10)*exp(-0.3*4^2); 若要显示变量 y 的值, 直接键入 y 即可: y y = -0.0045 在上例中, sin 是正弦函数, exp 是指数函数, 这些都是 MATLAB 常用到的数学函数. 下表即为 MATLAB 常用的基本数学函数及三角函数: 小整理 MATLAB 常用的基本数学函数 abs(x) 实数的绝对值或复数的模 angle(z) 复数 z 的相角(Phase angle) sqrt(x) 开平方 real(z) 复数 z 的实部 imag(z) 复数 z 的虚部 conj(z) 复数 z 的共轭复数 round(x) 四舍五入至最近整数 fix(x) 无论正负, 舍去小数至最近整数 floor(x) 大于 x 的最近整数 ceil(x) 小于 x 的最近整数 rat(x) 将实数 x 化为分数表示 rats(x) 将实数 x 化为多项分数展开 sign(x) 符号函数 (Signum function) 当x<0 时, sign(x)=-1; 当x=0 时, sign(x)=0; 当x>0 时, sign(x)=1. rem(x,y) 求x除以 y 的余数 gcd(x,y) 整数 x 和y的最大公因数 lcm(x,y) 整数 x 和y的最小公倍数 exp(x) 自然指数 ex 68 pow2(x) 2 的指数 2 x log(x) 以e为底的对数, 即自然对数或 ln(x) log2(x) 以2为底的对数 log2(x) log10(x) 以10 为底的对数 log10(x) 小整理 MATLAB 常用的三角函数 sin(x) 正弦函数 cos(x) 余弦函数 tan(x) 正切函数 asin(x) 反正弦函数 acos(x) 反余弦函数 atan(x) 反正切函数 atan2(x,y) 四象限的反正切函数 sinh(x) 双曲正弦函数 cosh(x) 双曲余弦函数 tanh(x) 双曲正切函数 asinh(x) 反双曲正弦函数 acosh(x) 反双曲余弦函数 atanh(x) 反双曲正切函数 变量也可用来存放向量或矩阵, 并进行各种运算, 如下例的行向量(Row vector)运算: x = [1 3 5 2]; y = 2*x+1 y = 3 7 11 5 小提示 变量命名的规则: 第一个字母必须是英文字母, 字母间不可留空格, 最多只能有 19 个字母, MATLAB 会忽略多余字母. 我们可以随意更改、增加或删除向量的元素: y(3) = 2 %更改第三个元素 y = 3 7 2 5 y(6) = 10 %加入第六个元素 y = 3 7 2 5 0 10 y(4)删除第四个元素, y = 3 7 2 0 10 在上例中, MATLAB 会忽略所有在百分比符号(%)之后的文字, 视为程序的注解(Comments). MATLAB 亦可取出向量的一个元素或一部份来做运算: x(2)*3+y(4) %取出 x 的第二个元素和 y 的第四个元素来做运算 ans = 9 y(2:4)-1 %取出 y 的第二至第四个元素来做运算 ans = 6 1 -1 69 在上例中, 2:4 代表一个由 2、3、4 组成的向量, 同样的方法可用于产生公差为 1 的等差数列: x = 7:16 x = 7 8 9 10 11 12 13 14 15 16 若不希望公差为 1, 则可将所需公差直接置于 4 与13 之间: x = 7:3:16 %公差为 3 的等差数列 x = 7 10 13 16 事实上, 我们可利用 linspace 来产生任意的等差数列: x = linspace(4, 10, 6) %等差数列: 首项为 4,末项为 10,项数为 6 x = 4.0000 5.2000 6.4000 5.6000 8.8000 10.0000 若对 MATLAB 函数用法有疑问, 可随时使用 help 来寻求在线帮助(on-line help). 将列向量转置(Transpose)后, 即可得到行向量(Column vector): z = x' z = 4.0000 5.2000 6.4000 5.6000 8.8000 10.0000 不论是行向量或列向量, 我们均可用相同的函数找出其元素个数、最大值、最小值等: length(z) %z 的元素个数 ans = 6 max(z) %z 的最大值 ans = 10 min(z) %z 的最小值 ans = 4 小整理 适用于向量的常用函数有 min(x) 向量 x 的元素的最小值 max(x) 向量 x 的元素的最大值 mean(x) 向量 x 的元素的平均值 median(x) 向量 x 的元素的中位数 std(x) 向量 x 的元素的标准差 diff(x) 向量 x 的相邻元素的差 sort(x) 对向量 x 的元素进行排序(Sorting) length(x) 向量 x 的元素个数 norm(x) 向量 x 的欧氏(Euclidean)长度 sum(x) 向量 x 的元素总和 prod(x) 向量 x 的元素总乘积 70 cumsum(x) 向量 x 的累计元素总和 cumprod(x) 向量 x 的累计元素总乘积 dot(x, y) 向量 x 和y的内积 cross(x, y) 向量 x 和y的外积 (大部份的向量函数也可适用于矩阵) 若要输入矩阵, 则必须在每一列结尾加上分号(;), 如下例: A = [1 2 3 4; 5 6 7 8; 9 10 11 12]; A 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 = 1 2 3 4 5 5 6 5 8 6 9 10 11 12 5 A(:, 2)删除第二行(: 代表所有列) A = 1 3 4 5 5 5 8 6 9 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 这几种矩阵处理的方式可以相互迭代运用, 产生各种意想不到的效果, 就看各位的巧思和创意. 小提示 在MATLAB 的内部资料结构中, 每一个矩阵都是一个以行为主(Column-oriented)的阵列(Array) 71 因此对于矩阵元素的存取, 我们可用一维或二维的索引(Index)来定址. 举例来说, 在上述矩阵 A 中, 位于第二列、第三行的元素可写为 A(2,3) (二维索引)或A(6)(一维索引, 即将所有直行进行堆叠后的 第六个元素). 此外, 若要重新安排矩阵的形状, 可用 reshape 命令: B = reshape(A, 4, 2) %4 是新矩阵的列数, 2 是新矩阵的行数 B = 5 8 9 12 5 6 11 5 小提示 A(:)就是将矩阵 A 每一列堆叠起来, 成为一个行向量, 而这也是 MATLAB 变量的内部储存方式. 以前例而言, reshape(A, 8, 1)和A(:)同样都会产生一个 8*1 的矩阵. . MATLAB 可在同时执行数个命令, 只要以逗号或分号将命令隔开: x = sin(pi/3); y = x^2; z = y*10, z = 5.5000 若一个数学运算是太长, 可用三个句点将其延伸到下一行: z = 10*sin(pi/3)* ... sin(pi/3); MATLAB 有些永久常数(Permanent constants), 使用者可直接取用, 例如: pi ans = 3.1416 下表即为 MATLAB 常用到的永久常数. 小整理 MATLAB 的永久常数 i 或j基本虚数单位(即1?)eps 系统的浮点(Floating-point)精确度 inf 无限大, 例如 1/0 nan 或NaN 非数值(Not a number), 例如 0/0 pi 圆周率(= 3.1415926...) realmax 系统所能表示的最大数值 realmin 系统所能表示的最小数值 nargin 函数的输入参数个数 nargin 函数的输出参数个数 §1.3 基本平面绘图命令 MATLAB 不但擅长于矩阵相关的数值运算, 也适合用在各种科学可视化表示(Scientific visualization). 本节将介绍 MATLAB 基本 xy 平面及 xyz 空间的各项绘图命令, 包含一维曲线及二维 曲面的绘制、输出及存档. plot 是绘制一维曲线的基本函数, 但在使用此函数之前, 我们需先定义曲线上每一点的 x 及y座72 标. 下例可画出一条正弦曲线: close all; x=linspace(0, 2*pi, 100); %100 个点的 x 座标 y=sin(x); %对应的 y 座标 plot(x,y); 小整理: MATLAB 基本绘图函数 plot x 轴和 y 轴均为线性刻度(Linear scale) loglog x 轴和 y 轴均为对数刻度(Logarithmic scale) semilogx x 轴为对数刻度, y 轴为线性刻度 semilogy x 轴为线性刻度, y 轴为对数刻度 若要画出多条曲线, 只需将坐标对依次放入 plot 函数即可: plot(x, sin(x), x, cos(x)); 若要改变颜色, 在座标对后面加上相关字串即可: plot(x, sin(x), 'c', x, cos(x), 'g'); 若要同时改变颜色及图线型态(Line style), 也是在座标对后面加上相关字串即可: plot(x, sin(x), 'co', x, cos(x), 'g*'); 小整理 plot 绘图函数的参数 字元 颜色 字元 图线型态 y 黄色 . 点k黑色 o 圆w白色 x x b 蓝色 + + g 绿色 * * r 红色 - 实线 c 亮青色 : 点线 m 锰紫色 -. 点虚线 -- 虚线 图形完成后, 我们可用 axis([xmin,xmax,ymin,ymax])函数来调整图轴的范围: axis([0, 6, -1.2, 1.2]); 此外, MATLAB 也可对图形加上各种注解与处理: xlabel('Input Value'); %x 轴注解 73 ylabel('Function Value'); %y 轴注解 title('Two Trigonometric Functions'); %图形标题 legend('y = sin(x)','y = cos(x)'); %图形注解 text(x,y, 'str'); %在坐标(x,y)处显示字符串 str grid on; %显示格线 我们可用 subplot 来同时画出数个小图形于同一个视窗之中: subplot(2,2,1); plot(x, sin(x)); subplot(2,2,2); plot(x, cos(x)); subplot(2,2,3); plot(x, sinh(x)); subplot(2,2,4); plot(x, cosh(x)); MATLAB 还有其他各种二维绘图函数, 以适合不同的应用, 详见下表. 小整理: 其他各种二维绘图函数 bar 长条图 errorbar 图形加上误差范围 fplot 较精确的函数图形 polar 极座标图 hist 累计图 rose 极座标累计图 stairs 阶梯图 stem 针状图 fill 实心图 feather 羽毛图 compass 罗盘图 quiver 向量场图 §1.4 基本三维绘图命令 plot3 命令将绘制二维图形的函数 plot 的特性扩展到三维空间. 函数格式除了包括第三维的信息 (比如 Z 方向)之外, 与二维函数 plot 相同. plot3 一般语法调用格式是 plot3(x1,y1,z1,S1,x2,y2,z2,S2,…) 这里 xn,yn 和zn 是向量或矩阵, Sn 是可选的字符串, 用来指定颜色、标记符号和线形. 总的来说, plot3 可用来画一个单变量的三维函数. 如下为一个三维螺旋线例子: t=0:pi/50:10*pi; plot3(sin(t),cos(t),t) title('Helix'),xlabel('sint(t)'),ylabel('cos(t)'),zlabel('t') text(0,0,0, 'Origin') grid v = axis v = -1 1 -1 1 0 40 74 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 0 10 20 30 40 sint(t) Origin Helix cos(t) t 从上例可明显看出, 二维图形的所有基本特性在三维中仍都存在. axis 命令扩展到三维, 只是返 回Z轴界限(0 和40), 在数轴向量中增加两个元素. 函数 zlabel 用来指定 z 轴的数据名称, 函数 grid 在图底绘制三维网格. 函数 test(x,y,z, 'string')在由三维坐标 x,y,z 所指定的位置放一个字符串. 另外, 子图和多图形窗口可以直接应用到三维图形中. 通过指定 plot3 命令的多个参数或使用 hold 命令, 可以把多条直线或曲线重叠画出. 例如, 增加 维数的 plot3 命令可以使多个二维图形沿一个轴排列起来, 而不是直接将二维图形叠到另一个的上 面. x=linspace(0,3*pi); % x-axis data z1=sin(x); % plot in x-z plane z2=sin(2*x); z3=sin(3*x); y1=zeros(size(x)); % spread out along y-axis y3=zeros(size(x)); % by giving each diffent y-axis values y2=y3/2; plot3(x,y1,z1,x,y2,z2,x,y3,z3); grid,xlabel('x-axis'),ylabel('y-axis'),zlabel('z-axis') title('sin(x),sin(2x),sin(3x)') 0 2 4 6 8 10 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 x-axis sin(x),sin(2x),sin(3x) y-axis z-axis 上述图形也可以沿另外方向堆列. plot3(x,z1,y1,x,z2,y2,x,z3,y3) grid,xlabel('x-axis'),ylabel('y-axis'),zlabel('z-axis') 75 title('sin(x),sin(2x),sin(3x)') 0 2 4 6 8 10 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 x-axis y-axis z-axis sin(x),sin(2x),sin(3x) §1.5 流程控制 计算机编程语言和可编程计算器提供许多功能, 它允许你根据决策结构控制命令执行流程. 如 果你以前已经使用过这些功能, 对此就会很熟悉. 相反, 如果不熟悉控制流, 本章材料初看起来或许 复杂些. 如果这样, 就慢慢来. 控制流极其重要, 因为它使过去的计算影响将来的运算. MATLAB 提供三种决策或控制流结构. 它们是: for 循环, while 循环和 if-else-end 结构. 由于这些结构经常包含大量的 MATLAB 命令, 故经 常出现在 M 文件中, 而不是直接加在 MATLAB 提示符下. 5.1 for 循环 for 循环允许一组命令以固定的和预定的次数重复. for 循环的一般形式是: for x = array {commands} end 在for 和end 语句之间的{commands}按数组中的每一列执行一次. 在每一次迭代中, x 被指定为 数组的下一列, 即在第 n 次循环中, x=array(:, n). 例如, for n=1:10 x(n)=sin(n*pi/10); end x 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 换句话, 第一语句是说: 对n等于 1 到10, 求所有语句的值, 直至下一个 end 语句. 第一次通过 for 循环 n=1, 第二次 n=2, 如此继续, 直至 n=10. 在n=10 以后, for 循环结束, 然后求 end 语句后面的 任何命令值, 在这种情况下显示所计算的 x 的元素. for 循环的其它重要方面是: 1. for 循环不能用 for 循环内重新赋值循环变量 n 来终止. for n=1:10 76 x(n)=sin(n*pi/10); n=10; end x 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 2. 语句 1:10 是一个标准的 MATLAB 数组创建语句. 在for 循环内接受任何有效的 MATLAB 数组. data=[3 9 45 6; 7 16 -1 5] data = 3 9 45 6 7 16 -1 5 for n=data x=n(1)-n(2) end x = -4 x = -7 x = 46 x = 1 3. for 循环可按需要嵌套. for n=1:5 for m=5:-1:1 A(n,m)=n^2+m^2; end disp(n) end 1 2 3 4 5 A A = 2 5 10 17 26 5 8 13 20 29 10 13 18 25 34 17 20 25 32 41 77 26 29 34 41 50 4. 当有一个等效的数组方法来解给定的问题时, 应避免用 for 循环. 例如, 上面的第一个例子可 被重写为 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 两种方法得出同样的结果, 而后者执行更快, 更直观, 要求较少的输入. 5. 为了得到最大的速度, 在for 循环(while 循环)被执行之前, 应预先分配数组. 例如, 前面所考 虑的第一种情况, 在for 循环内每执行一次命令, 变量 x 的大小增加 1. 迫使 MATLAB 每通过一次循 环要花费时间对 x 分配更多的内存. 为了消去这个步骤, for 循环的例子应重写为 x=zeros(1,10); % preallocated memory for x for n=1:10 x(n)=sin(n*pi/10); end 现在, 只有 x(n)的值需要改变. 5.2 while 循环 与for循环以固定次数求一组命令的值相反, while循环以不定的次数求一组语句的值. while循环 的一般形式是: while expression {commands} end 只要在表达式里的所有元素为真, 就执行 while 和end 语句之间的{commands}. 通常, 表达式的求值 给出一个标量值, 但数组值也同样有效. 在数组情况下, 所得到数组的所有元素必须都为真. 考虑下 列例子: num=0;EPS=1; while (1+EPS)>1 EPS=EPS/2; num=num+1; end num num = 53 EPS=2*EPS EPS = 2.2204e-016 这个例子表明了计算特殊 MATLAB 值EPS 的一种方法, 它是一个可加到 1, 而使结果以有限精 度大于 1 的最小数值. 这里我们用大写 EPS, 因此 MATLAB 的EPS 的值不会被覆盖掉. 在这个例子 里, EPS 以1开始. 只要(1+EPS)>1 为真(非零), 就一直求 while 循环内的命令值. 由于 EPS 不断地被 2 除, EPS 逐渐变小以致于 EPS+1 不大于 1. (记住, 发生这种情况是因为计算机使用固定数的数值来 表示数. MATLAB 用16 位, 因此, 我们只能期望 EPS 接近 10-16. ) 在这一点上, (1+EPS)>1 是假(零), 78 于是 while 循环结束. 最后, EPS 与2相乘, 因为最后除 2 使EPS 太小. 5.3 if-else-end 结构 很多情况下, 命令的序列必须根据关系的检验有条件地执行. 在编程语言里, 这种逻辑由某种 if-else-end 结构来提供. 最简单的 if-else-end 结构是: if expression {commands} end 如果在表达式中的所有元素为真(非零), 那么就执行 if 和end 语言之间的{commands}. 在表达式 包含有几个逻辑子表达式时, 即使前一个子表达式决定了表达式的最后逻辑状态, 仍要计算所有的 子表达式. 例如, apples=10; % number of apples cost=apples*25 % cost of apples cost = 250 if apples>5 % give 20% discount for larger purchases cost=(1-20/100)*cost; end cost cost = 200 假如有两个选择, if-else-end 结构是: if expression commands evaluated if True else commands evaluated if False end 在这里, 如果表达式为真, 则执行第一组命令; 如果表达式是假, 则执行第二组命令. 当有三个或更多的选择时, if-else-end 结构采用形式 if expression1 commands evaluated if expression1 is True elseif expression2 commands evaluated if expression2 is True elseif expression3 commands evaluated if expression3 is True elseif expression4 commands evaluated if expression4 is True elseif …… …… else commands evaluated if no other expression is True end 最后的这种形式, 只和所碰到的、 与第一个真值表达式相关的命令被执行; 接下来的关系表达式 不检验, 跳过其余的 if-else-end 结构. 而且, 最后的 else 命令可有可无. 现在我们知道了如何用 if-else-end 结构来决策, 就有可能提出一种合理的方法来跳出或中断 for 79 循环和 while 循环. EPS=1; for num=1:1000 EPS=EPS/2; if (1+EPS)<=1 EPS=EPS*2 break end end EPS = 2.2204e-016 num num = 53 这个例子演示了估算 EPS 的另一种方法. 在这种情况下, for 循环构造成要执行足够多的次数. if-else-end 结构检验要看 EPS 是否变得足够小. 如果是, EPS 乘2, break 命令强迫 for 循环提早结束, num=53. 在这个例子里, 当执行 break 语句时, MATLAB 跳到循环外下一个语句. 在现在情况下, 它返回 到MATLAB 的提示符并显示 EPS. 如果一个 break 语句出现在一个嵌套的 for 循环或 while 循环结构 里, 那么 MATLAB 只跳出 break 所在的那个循环, 不跳出整个嵌套结构. 5.4 小结 MATLAB 控制流功能可以概括如下: 控制流结构 for x = array commands end for 循环, 每一次迭代将 x 赋以数组的 第i列, 执行命令 while expression commands end while 循环, 只要表达式的所有元素 为真或非零, 执行命令. if expression commands end 简单的 if-else-end 结构, 若在表达式 中的所有元素是真值或非零, 执行 命令. if expression commands evaluated if True else expression commands evaluated if False end 具有两条路径的 if-else-end 结构, 若表达式为真或非零, 则执行一组 命令. 若表达式为假或零, 则执行另一组命令. if expression1 commands evaluated if expression1 is True elseif expression2 commands evaluated if expression2 is True if expression3 commands evaluated if expression3 is True elseif …… 最一般的 if-else-end 结构. 只执行与第一个真值表达式相关 的命令. 80 …… else commands evaluated if no other expression is True end break 结束 for 循环和 while 循环的执行 §1.6 M 文件 Matlab 允许将一系列将要运行的语句事先存放在后缀名为.m 的文件中, 然后在 Matlab 命令行中 输入该文件的主文件名即可. 当然, 要让 Matlab 能过文件名找得到文件并执行相应的语句, 必须把 文件放到 Matlab 的Current Directory 中, 该目录可以在 Matlab 界面中自己设定, 默认是安装目录下 的work 目录. M 文件是普通的文本文件, 可以用记事本或其它的文本编辑工具编写. Matlab 提供了 M 文件编辑工具, 可通过菜单 File\new\M-file 新建一个 M 文件, 也可以用 File\Open 打开一个已有的 M 文件进行修改. 编写 M 文件有两个目的. 一是将多个 Matlab 语句集成在一个文件, 然后通过一个命令一次运行. 这相当于批处理文件. 例如用记事本输入以下内容, 并以 test.m 为文件名保存到 Matlab 的Current Directory 中: n=50; b=3; a=0; x=linspace(a,b,n); e1=exp(-x.^2); e2=(x.^2).*exp(-x.^2); e3=x.*exp(-x.^2); e4=exp(-x); plot(x,e1,x,e2,x,e3,x,e4) 然后在 Matlab 命令行中 test, 可得到如下图形: 0 0.5 1 1.5 2 2.5 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 二是编写自定义函数. Matlab 的自定义函数以 M 文件形式存储在计算机中, 一个函数用一个 M 文件存储, 这类 M 文件称为函数 M 文件. 函数 M 文件的一般格式: function 函数值=函数名(参数) 81 函数体 function 是关键字, 表明下面定义一个函数. 函数值用来取得函数的返回值, 通常是在函数体中 对其进行赋值, 一个函数可以有多个函数值, 多个函数值用逗号隔开, 有时候也称函数值为输出参 数. 函数名必须与 M 文件的主文件名相同, 是Matlab 调用该函数的依据, 函数名最多可以有 31 个字 符, 并且与操作系统有关, 有些操作系统允许的字符数可能更少, 其命名规则与变量命名规则相同. 参数用来从函数外部获取数据, 可以有 0 个或多个参数, 多个参数之间用逗号隔开, 有时候也称参数 为输入参数. 函数体是实现该函数的功能的语句, 当最后一行被执行完毕或遇到return 语句时退出函 数, 并返回到 Matlab 命令行. 例如, 在记事本中输入以下内容, 并以文件名 sqpulse.m 存到 Matlab 的Current Directory 指定的目录中: function f=sqpulse(n,x) % 递归函数用于求和 % 1/2+2/pi cos(x pi)2sin(n pi/2) cos(n x pi) % 当n-->inf 这将等于阶跃的平方 if(n==1) f=1/2+2/pi*cos(x*pi); else f=2*sin(n*pi/2)/n/pi*cos(n*x*pi)+sqpulse(n-1,x); end 然后在 Matlab 命令行中输入 n=100; x=-3:0.1:3; plot(x,sqpulse(n,x)) %画出函数 sqpulse 的图形 可得到如下图形: -3 -2 -1 0 1 2 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 §1.7 曲线拟合与插值 在大量的应用领域中, 人们经常面临用一个解析函数描述数据(通常是测量值)的任务. 对这个问 题有两种方法. 在插值法里, 数据假定是正确的, 要求以某种方法描述数据点之间所发生的情况. 这 种方法在下一节讨论. 这里讨论的方法是曲线拟合或回归. 人们设法找出某条光滑曲线, 它最佳地 拟合数据, 但不必要经过任何数据点. 图6.1 说明了这两种方法. 标有'o'的是数据点; 连接数据点的 实线描绘了线性内插, 虚线是数据的最佳拟合. 82 6.1 曲线拟合 曲线拟合涉及回答两个基本问题: 最佳拟合意味着什么?应该用什么样的曲线?可用许多不同 的方法定义最佳拟合, 并存在无穷数目的曲线. 所以, 从这里开始, 我们走向何方?正如它证实的那 样, 当最佳拟合被解释为在数据点的最小误差平方和, 且所用的曲线限定为多项式时, 那么曲线拟 合是相当简捷的. 数学上, 称为多项式的最小二乘曲线拟合. 如果这种描述使你混淆, 再研究图 6.1. 虚线和标志的数据点之间的垂直距离是在该点的误差. 对各数据点距离求平方, 并把平方距离全加 起来, 就是误差平方和. 这条虚线是使误差平方和尽可能小的曲线, 即是最佳拟合. 最小二乘这个术 语仅仅是使误差平方和最小的省略说法. 0 0.2 0.4 0.6 0.8 1 -2 0 2 4 6 8 10 12 x y=f(x) Second Order Curve Fitting 图6.1 2 阶曲线拟合 在MATLAB 中, 函数 polyfit 求解最小二乘曲线拟合问题. 为了阐述这个函数的用法, 让我们以 上面图 6.1 中的数据开始. x=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1]; y=[-.447 1.978 3.28 6.16 7.08 5.34 5.66 9.56 9.48 9.30 11.2]; 为了用 polyfit, 我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶次或度. 如果我们选择 n=1 作为阶次, 得到最简单的线性近似. 通常称为线性回归. 相反, 如果我们选择 n=2 作为阶次, 得到一个 2 阶多项式. 现在, 我们选择一个 2 阶多项式. n=2; % polynomial order p=polyfit(x, y, n) p = -9.8108 20.1293 -0.0317 polyfit 的输出是一个多项式系数的行向量. 其解是 y = -9.8108x2 +20.1293x-0.0317. 为了 将曲线拟合解与数据点比较, 让我们把二者都绘成图. xi=linspace(0, 1, 100); % x-axis data for plotting z=polyval(p, xi); 为了计算在 xi 数据点的多项式值, 调用 MATLAB 的函数 polyval. plot(x, y, ' o ' , x, y, xi, z, ':' ) 画出了原始数据 x 和y, 用'o'标出该数据点, 在数据点之间, 再用直线重画原始数据, 并用点':'线, 画出多项式数据 xi 和z. xlabel(' x '), ylabel(' y=f(x) '), title(' Second Order Curve Fitting ') 83 将图作标志. 这些步骤的结果表示于前面的图 6.1 中. 多项式阶次的选择是有点任意的. 两点决定一直线或一阶多项式. 三点决定一个平方或 2 阶多 项式. 按此进行, n+1 数据点唯一地确定 n 阶多项式. 于是, 在上面的情况下, 有11 个数据点, 我们可 选一个高达 10 阶的多项式. 然而, 高阶多项式给出很差的数值特性, 人们不应选择比所需的阶次高 的多项式. 此外, 随着多项式阶次的提高, 近似变得不够光滑, 因为较高阶次多项式在变零前, 可多 次求导. 例如, 选一个 10 阶多项式 pp=polyfit(x, y, 10) ; format short e % change display format pp.' % display polynomial coefficients as a column ans = -4.6436e+005 2.2965e+006 -4.8773e+006 5.8233e+006 -4.2948e+006 2.0211e+006 -6.0322e+005 1.0896e+005 -1.0626e+004 4.3599e+002 -4.4700e-001 要注意在现在情况下, 多项式系数的规模与前面的 2 阶拟合的比较. 还要注意在最小 (-4.4700e-001)和最大(5.8233e+006)系数之间有 7 个数量级的幅度差. 将这个解作图, 并把此图与原 始数据及 2 阶曲线拟合相比较, 结果如何呢? zz=polyval(pp, xi); % evaluate 10th order polynomial plot(x, y, ' o ' , xi, z,xi, zz) % plot data xlabel(' x '), ylabel(' y=f(x) '), title(' 2nd and 10th Order curve Fitting ') 在下面的图 6.2 中, 原始数据标以'o', 2 阶曲线拟合是虚线, 10 阶拟合是实线. 注意, 在10 阶拟合 中, 在左边和右边的极值处, 数据点之间出现大的纹波. 当企图进行高阶曲线拟合时, 这种纹波现象 经常发生. 根据图 6.2, 显然, "越多就越好"的观念在这里不适用. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -2 0 2 4 6 8 10 12 14 16 18 x y=f(x) 2nd and 10th Order curve Fitting 图6.2 2 阶和 10 阶曲线拟合 84 6.2 一维插值 正如在前一节对曲线拟合所描述的那样, 插值定义为对数据点之间函数的估值方法, 这些数据 点是由某些集合给定. 当人们不能很快地求出所需中间点的函数值时, 插值是一个有价值的工具. 例如, 当数据点是某些实验测量的结果或是过长的计算过程时, 就有这种情况. 或许最简单插值的例子是 MATLAB 的作图. 按缺省, MATLAB 用直线连接所用的数据点以作图. 这个线性插值猜测中间值落在数据点之间的直线上. 当然, 当数据点个数的增加和它们之间距离的 减小时, 线性插值就更精确. 例如, x1=linspace(0, 2*pi, 60); x2=linspace(0, 2*pi, 6); plot(x1, sin(x1), x2, sin(x2) xlabel(' x '), ylabel(' sin(x) '), title(' Linear Interpolation ') 0 1 2 3 4 5 6 7 -1 -0.5 0 0.5 1 x sin(x) Linear Interpolation 图6.3 线性插值 图6.3 是sine 函数的两个图, 一个在数据点之间用 60 个点, 它比另一个只用 6 个点更光滑和更 精确. 如曲线拟合一样, 插值要作决策. 根据所作的假设, 有多种插值. 而且, 可以在一维以上空间中 进行插值. 即如果有反映两个变量函数的插值, z=f(x, y), 那么就可在 x 之间和在 y 之间, 找出 z 的中 间值进行插值. MATLAB 在一维函数 interp1 和在二维函数 interp2 中, 提供了许多的插值选择. 其中 的每个函数将在下面阐述. 为了说明一维插值, 考虑下列问题, 12 小时内, 一小时测量一次室外温度. 数据存储在两个 MATLAB 变量中. hours=1:12; % index for hour data was recorded temps=[5 8 9 15 25 29 31 30 22 25 27 24]; % recorded temperatures plot(hours, temps, hours, temps,view temperatures title(' Temperature ') xlabel(' Hour '), ylabel(' Degrees Celsius ') 85 0 2 4 6 8 10 12 5 10 15 20 25 30 35 Hour Degrees Celsius Temperature 图6.4 在线性插值下室外温度曲线 正如图 6.4 看到的, MATLAB 画出了数据点线性插值的直线. 为了计算在任意给定时间的温度, 人们可试着对可视的图作解释. 另外一种方法, 可用函数 interp1. t=interp1(hours, temps, 9.3) % estimate temperature at hour=9.3 t = 22.9000 t=interp1(hours, temps, 4.7) % estimate temperature at hour=4.7 t = 22 t=interp1(hours, temps, [3.2 6.5 5.1 11.7]) % find temp at many points! t = 10.2000 30.0000 30.9000 24.9000 interp1 的缺省用法是由 interp1(x, y, xo)来描述, 这里 x 是独立变量(横坐标), y 是应变量(纵坐标), xo 是进行插值的一个数值数组. 另外, 该缺省的使用假定为线性插值. 若不采用直线连接数据点, 我们可采用某些更光滑的曲线来拟合数据点. 最常用的方法是用一 个3阶多项式, 即3次多项式, 来对相继数据点之间的各段建模, 每个 3 次多项式的头两个导数与该 数据点相一致. 这种类型的插值被称为 3 次样条或简称为样条. 函数 interp1 也能执行 3 次样条插值. t=interp1(hours, temps, 9.3, ' spline ') % estimate temperature at hour=9.3 t = 21.8577 t=interp1(hours, temps, 4.7, ' spline ') % estimate temperature at hour=4.7 t = 22.3143 t=interp1(hours, temps, [3.2 6.5 5.1 11.7], ' spline ') t = 9.6734 86 30.0427 31.1755 25.3820 注意, 样条插值得到的结果, 与上面所示的线性插值的结果不同. 因为插值是一个估计或猜测 的过程, 其意义在于, 应用不同的估计规则导致不同的结果. 一个最常用的样条插值是对数据平滑. 例如, h=1:0.1:12; % estimate temperature every 1/10 hour t=interp1(hours, temps, h, ' spline ') ; plot(hours, temps,hours, temps,h, t) % plot comparative results title(' Springfield Temperature ') xlabel(' Hour '), ylabel(' Degrees Celsius ') 在图 6.5 中, 虚线是线性插值, 实线是平滑的样条插值, 标有' + '的是原始数据. 如要求在时间轴 上有更细的分辨率, 并使用样条插值, 我们有一个更平滑、但不一定更精确地对温度的估计. 尤其应 注意, 在数据点, 样条解的斜率不突然改变. 作为这个平滑插值的回报, 3 次样条插值要求更大量的 计算, 因为必须找到 3 次多项式以描述给定数据之间的特征. 0 2 4 6 8 10 12 5 10 15 20 25 30 35 Hour Degrees Celsius Springfield Temperature 图6.5 在不同插值下室外温度曲线 6.3 小结 下面的表 6.1 总结了在 MATLAB 中所具有的曲线拟合和插值函数. 表6.1 曲线拟合和插值函数polyfit(x, y, n) 对描述 n 阶多项式 y=f(x)的数据, 进行最小二乘曲线拟合 interp1(x, y, xo) 1 维线性插值 interp1(x, y, xo, ' spline ') 1 维3次样条插值 interp1(x, y, xo, ' cubic ') 1 维3次插值 interp2(x, y, Z, xi, yi) 2 维线性插值 interp2(x, y, Z, xi, yi, ' cubic ') 2 维3次插值 interp2(x, y, Z, xi, yi, ' nearest ') 2 维最近邻插值 87 §1.8 求解非线性规划 非线性规划的标准型为: min F(X) s.t.G(X) ≤0 VLB≤X≤VUB 其中 X∈En, G(X)∈Em . 用Matlab 求解, 分两步: 1.首先建立函数形式的 M 文件 fun.m: function[f,g]=fun(X); f=F(X); g=[G1(X); …; Gm(X)] 2.命令的基本格式有三个: (1)X=constr('fun', x0 ) (2)X=constr('fun', x0, options ) (3)X=constr('fun', x0, options, VLB, VUB ) 其中 x0, X 为迭代初值, options 是参数说明向量. 例 求解以下模型 1 2 2 1 2 1 2 2 1 2 1 2 1 2 1 2 min ( ) (4 2 4 2 1) 0 . . 1.5 0 10 0 x f x e x x x x x x x s t x x x x x x + = ? ? + ? ? ≤ ? ? ? ? ≤ ? 解: (1)先建立 M-文件 fun1.m: function[f,g]=fun(X); f=exp(x(1))﹡(4﹡x(1)^2+2﹡x(2)^2+4﹡x(1)﹡x(2)+2﹡x(2)+1); g=[x(1)+x(2); 1.5+x(1)﹡x(2)-x(1)-x(2); -x(1)﹡x(2)-10]; (2)再输入以下命令: x0=[-2,2]; options(13)=1; (说明 g 中第一个约束为等式约束) [x, options]= constr('fun1', x0, options ) options(8) (表示输出在最后极值点的函数值) (3)运算结果: x= -1.2247 1.2247 ans= 1.8951 第二节 Maple 现在的计算机不仅可以进行数值计算, 还可以进行符号运算. 数值计算就是处理的数据和处理 后的结果都是数值, 而符号运算则是处理的数据和处理后的结果都是符号. 所谓符号, 可以是字母、 公式, 也可以是数字. 符号计算实际上就是数学演算、数学推理和数学证明. 专门用来完成这些演 88 算、推理和证明的软件, 通常称为数学软件. 当今比较成熟的数学软件有 MatLab、Maple 和mathematic. MatLab 主要是进行矩阵的运算, 具有很强的数值运算能力; 而Maple 和mathematic 则擅 长符号运算, 例如公式的推导、变形、方程的符号解、微积分运算等等; 它们各有自己的优点, 并列 为三大数学软件. §1.1 初识 Maple 1. Maple 简介 Maple 是一个优秀的国际流行的数学软件系统, 它提供了广泛的数学计算功能, 主要包括三个 方面: 符号演算、二维和三维图形描绘以及数值计算. Maple 能解决如下问题: 初等数学中的有理数运算(精确计算和大有理数计算), 代数式运算, 初等函数变换, 映射与集合 运算, 方程求根等. 复数运算, 向量代数, 平面和空间解析几何, 线性代数中矩阵的各种运算和变换, 线性方程组 求解. 微积分中的函数的极限, 导数, 积分, 泰勒级数、傅立叶级数展开, 以及数列与级数求和, 求积 等. 微分方程的解析解, 数值解以及通过方向场图形得出近似解. 此外, Maple 还提供了解决诸如图论、数论、群论、组合数学、概率统计、特殊函数(运算微积) 和其他应用数学(如金融数学)方面的问题的软件包. 总之, 凡是能用公式解决的问题, Maple 都能解 决; 它能有效地辅助人们在数学计算方面的繁琐事务, 进而去从事更加智能化和深层次的问题的研 究. Maple 是加拿大 Waterloo 大学符号计算小组于 1980 年开始研制, 经过多年的应用和开发, 至今 发展到 Maple9. 由于枫树是加拿大的国树, Maple 的开发者为了表示自己的爱国之情, 所以把枫叶作 为Maple 的图标和标志. Maple 系统包含有三部分: 用户界面 Iris, 基本部分(内核)Kernel 和外部库 Liberary: Iris. 语法分析(对用户的输入进行分析), 主要功能: 表达式显示, 图形处理, 特殊用户界面. Kernel. 翻译器(将用户输入翻译成机器代码), 主要功能: 内存管理, 整数域, 有理数域, 实数域, 扩展有理数域的计算和处理. Liberary. 外部库函数, 包括应用软件包和在线帮助. Maple 操作界面 用户界面 Iris 和基本部分 Kernel 是Maple 的主体, 由C++语言编写而成, 用户启用使用 Maple 时就将它们调入内存, 直接可以执行. 外部库 Liberary 的库函数由 Maple 语言编写而成, 在用户需要 时才调入内存. Maple 还提供了方便的帮助系统, 只要键入相关的关键字, 按下 F1 键, 或者在工作窗口中输入"? 89 关键字"再按 Enter, 就可以得到详细的说明和例子; 通常, 每一个在线帮助系统的画面的最下方都有 有一些链接(links), 提供了一些其它相关的主题. 如果要查询哪个主题, 按一下该链接即可显示相关 性的主题的画面. 2. 第一次使用 Maple 启动 Maple 之后, Maple 会自动开启一个新的"工作窗口". 每一个工作窗口即为一个文件, 而Maple 用一种特殊的文件格式来存储工作窗口中的内容, 其文件名后缀为.mws (Maple WorkSheet). 假设要求解方程式 x2 +ax-b=0 的根, 可以在工作窗口中键入: solve(x^2+a*x-b=0,{x}); 再按 Enter 键来执行这个简单的运算, 此时您会得到如下的画面: 注意事项: 1. 上图中的提示符号(>)为Maple 自动产生. 当您输入算式时, 一定不要把提示符号也敲进去, 否则会有错误的消息产生. 2. Maple 使用":="作为赋值号而不是"=". 3. Maple 表达式的结束符是";"和":", 以";"结束时 Maple 显示结果, 而以":"结束则不显示结果. 命令行最左边的"["所框的范围代表一个工作组, Maple 每收到一个执行指令(按Enter 键或点击 按钮)就执行当前光标在的工作内的所有命令. 也就是说, 在同一个工作组内的命令总是同时执 行的. 如果想一次执行一个文件内的所有命令, 可以点击 按钮. 通常 Maple 在执行完最后一条命令后在后面产生一个空的命令输入行. 如果想人为增加, 可用 鼠标点击工具栏上的 按钮. 4. 因为按 Enter 是执行命令, 想强迫换行必须按 Shift+Enter. 当执行的语句很多以后, 前面变量的值会影响到后面的语句而得到不是您所期望的结果, 可以 通过执行 restart;或点击 按钮来清除所有用户定义变量的值. §1.2 初等数学 1. 阶乘!和factorial 格式: n! 或factorial(n) ——求n的阶乘. 2. 绝对值 abs 格式: abs(x) ——求x的绝对值. 3. 质数测试 isprime 及求质数 nextprime 90 格式: isprime(n) ——测试整数 n 是否是质数. nextprime(n) ——求出比 n 大的最小质数. 4. 最大数 max 和最小数 min 格式: max(x1, x2,相当于 max{x1, x2, ...}, 即求 x1, x2, ...的最大数. min(x1, x2,求x1, x2, ...的最小数. 5. 随机数 rand 格式: rand()——返回 12 位的非负随机整数. rand(a..b) ——生成一个函数, 该函可以产生区间[a,b]内的随机整数. 6. 多项式插值函数 interp 格式: interp(x, y, v) ——函数 interp(x, y, v)根据下列独立点(x[1],y[1]), (x[2],y[2]x[n],y[n]) 计算小于或等于 n 次的插值多项式. 例with(plots): p:=interp([0,1,2,3],[0,3,1,3],x); poin:=pointplot({[0,0],[1,3],[2,1],[3,3]},color=green): #作数据点图, 只把数据用 poin 保存, 暂不 显示 line:=plot(p,x=-0.5..3.5): #作多项式 p 的图形, 暂不显示 display(poin,line); #显示图形 := p ? + 3 2 x3 7 x2 17 2 x 7. 排序 sort 格式: 1. sort(L) 若L是表, 则sort(L)将L的元素按升序排序; 若L是代数表达, 则sort(L)将L按降幂排列; 2. sort(L,F) F 是两个参数的布尔过程, 按F指定方法进行降幂排列. 8. 求和 sum、add 与求积 product、mul 格式: sum(f, k=m..n) 计算 ∑ = n m k f . add(f, k=m..n) 计算 ∑ = n m k f . product(f, k=m..n) 计算∏ = n m k f . 91 mul(f, k=m..n) 计算∏ = n m k f . 9. 几个用于表达式化简的函数——simplify、factor、combine、expand 格式: simplify(expr) 对表达式 expr 进行化简, 化简一词在数学上是含糊不清的, 但Maple 常常 能得到令人满意的结果. 如果用 simplify 化简仍达不到要求, 可以考虑结合其它几个化 简函数一起使用. factor(expr) 对表达式 expr 进行因式分解. combine(expr) 对表达式 expr 进行合并. expand(expr) 对表达式 expr 进行展开. 这些函数是非常有用的, 如果某个命令求得的结果非常的冗长, 您应该尝试用这些命令对它进 行化简. 这些命令的用法在此并没有完全列出, 详细请查阅 Maple 系统自带的帮助. 例simplify(cos(x)^5 + sin(x)^4 + 2*cos(x)^2 - 2*sin(x)^2 - cos(2*x)); + ( ) cos x 5 ( ) cos x 4 combine(exp(sin(a)*cos(b))*exp(cos(a)*sin(b))); e ( ) sin + a b factor(x^3+y^3); ( ) + x y ( ) ? + x2 x y y2 expand(sin(x+y)); + ( ) sin x ( ) cos y ( ) cos x ( ) sin y 10. 求表达式值 eval、evalf、evalm、evalb、value (1) 格式: eval(expr) 求出表达式 expr 的最终值. eval(expr, x=a) 求出表达式 expr 在x=a 处的值. (2) 格式: evalf(expr) 以默认的精度求出表达式 expr 的近似值. evalf[n](expr) 求出表达式 expr 近似值, 并且精确到小数点后面 n 位. (3) 格式: evalm(matrix expr) 求出矩阵 表达式 matrix expr 的值. Maple 为了提高运行速度, 有些变量并不求出最终的值, 而是以变量的形式参与, 矩阵运算尤其如此. 为了求出矩阵的值, 使用 evalm. (4) 格式: evalb(Boolean expr) 求出布尔表达式 matrix expr 的值 true 或false, 如果无法判断返 回原表达式. 布尔表达式由关系符>(大于)、 >=(大于等于)、 =(等于)、 <=(小于等于)、 <(小于)、 <>(不 等于)及逻辑操作符 and(与)、or(或)、not(非)将变量或表达式连接而成. (5) 格式: value(expr) 求出表达式 expr 的值. Maple 的极限、求导、积分、求和、求积等运算中, 都有一种命令是只写出求表达 式而没有真正计算, 这些命令通常以大写字母开头, 例如极限 Limit, 求导 Diff, 积分Int, 求和 Sum, 求积 Product, 如果要求出其值, 一是使用小写字母开头的命令, 再就是使用 value(). 92 §1.3 方程求解 1. 解方程 solve 格式: solve(eqns, vars) 说明: solve(eqns, vars)用来求解方程或方程组 eqns, vars 指出未知数. eqns 可以是线性方程、非线性方程, 可以是多项式、三解函数、指数和绝对值方程, 甚至可以 是不等式. 如果方程不存在显式解或显式解过于复杂, 则返回 RootOf 形式的解, RootOf(expr)的意思是 expr 的解. 解绝对值方程,函数 solve 假定 abs 的参数是实数值. 有时, Maple 不能找到解, 并不表明解不存在, 只是 Maple 不能找到它. 这时可以尝试用 fsolve 求它的近似解. 例eqns := {u+v+w=1, 3*u+v=3, u-2*v-w=0}; #解三元线性方程组 := eqns { } , , = + + u v w 1 = + 3 u v 3 = ? ? u 2 v w 0 sols := solve( eqns ); := sols { } , , = w -2 5 = u 4 5 = v 3 5 eq := x^4-5*x^2+6*x=2; #求解高次多项式的根 := eq = ? + x4 5 x2 6 x 2 solve(eq,x); , , , 1 1 ? + 1 3 ? ? 1 3 solve(x-sin(x), x); #求解三角函数议程 0 sovle(x^2+x>5, x); #求解不等式 ( ) sovle , < 5 + x2 x x solve(abs(x)>=5, x); #求解绝对值不等式 , ( ) RealRange , 5 ∞ ( ) RealRange , ?∞ -5 2. 求方程近似解 fsolve 格式: fsolve(eqns, vars) 说明: 对一般的方程, fsolve 计算出一个实根, 对多项式计算出所有实根(非复数). 有时用 solve 求方程的解的结果非常冗长或者是求不出(但解存在), 用fsolve 往往能得到满意的 结果. 例solve( tan(sin(x))=1, x ); ? ? ? ? ? ? ? ? arcsin 1 4 π fsolve( tan(sin(x))=1, x ); .9033391108 poly := 23*x^5 + 105*x^4 - 10*x^2 + 17*x: solve( poly, x); 93 0 ( ) RootOf , + ? + 23 _Z4 105 _Z3 10 _Z 17 = index 1 , , ( ) RootOf , + ? + 23 _Z4 105 _Z3 10 _Z 17 = index 2 , ( ) RootOf , + ? + 23 _Z4 105 _Z3 10 _Z 17 = index 3 , ( ) RootOf , + ? + 23 _Z4 105 _Z3 10 _Z 17 = index 4 fsolve( poly, x, -1..1 ); -4.536168981, -.6371813185, 0. 3. 方程整数解 isolve 格式: isolve(eqns, vars) 说明: isolve(eqns, vars)求出方程 eqns 的整数解. 返回解中若含有变量_N1, _N2, 其含义是通用变量. 如果返回 NULL, 则说明 Maple 找不到解, 但并没有保证真的无解. 例isolve(3*x-4*y=7); { } , = x + 5 4 _Z1 = y + 2 3 _Z1 4. 解模方程 msovle 格式: msolve(eqns, m) msolve(eqns, vars, m) 说明: msolve(eqns, m)求方程的所有求知数关于模 m 的整数解. 解里面含有通用变量_NN1, ..., _NNn 表示整数参数. msolve(eqns, vars, m)求方程指定求知数关于模 m 的整数解. 例msolve({3*x-4*y=1,7*x+y=2},19); { } , = y 11 = x 15 msolve(x^2=3,5); msolve(2^i=3,19); { } = i + 13 18 _Z1~ msolve(8^j=2,x,17); { } = j + 3 8 x 5. 解递推方程 rsovle 格式: rsolve(eqns, fcns) 说明: rsolve(eqns, fcns)解出递推关系指定的方程 eqns 中的 fcns, 返回 fcns 的通项表达式. 例rsolve({y(n)=n+y(n-1), y(0)=0}, y(n)); ? ? ( ) + n 1 ? ? ? ? ? ? ? ? + 1 2 n 1 n 1 simplify(%); + 1 2 n2 1 2 n rsolve({f(n) = f(n-1) + f(n-2), f(0)=1, f(1)=1},{f(n)}); 94 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ( ) f n + 2 5 5 ? ? ? ? ? ? ? ? 2 1 ? + 1 5 n ? + 1 5 2 5 5 ? ? ? ? ? ? ? ? ?2 1 + 5 1 n + 5 1 §1.4 线性代数(矩阵运算) 1. 数组定义 array 格式: array(bounds, list) 说明: 参数中至少有范围 bounds 或初值表 list. Maple 数组的下标可以是负数和零. array(bounds)生成一个未指定元素值的数组, array(list)直接生成数组, 元素值由 list 指定. array 可以定义多维数组, 形如 A:=array(1..10, 1..10)定义二维数组 A10*10,多维数组类似. 一维数组的list 形如[a1,a2,...,an], 二维数组的list 形如[[a11,a12,...,a1n],[a21,a22,...,a2n], ..., [am1,am2,...,amn]], 以此类推. Maple 以索引来表示数组. Maple 用A[1,1]来取得二维数组 A 的元素, 而不是 A[1][1]. C 语言中限定同一个数组中的元素必须具有相同的数据类型, Maple 并没有此限制. 例A:= array(-1..1, -1..1); := A ( ) array , , .. -1 1 .. -1 1 [ ] A[-1,0]:=1; := A , -1 0 1 A[1,1]:=2*x^2+x; := A , 1 1 + 2 x2 x lprint(eval(A)); array(-1 .. 1, -1 .. 1,[(-1, 0)=1,(1, 1)=2*x^2+x]) B:= array(0..1,0..1,[[1,2],[x,x^2]]); B .. 0 1 .. 0 1 array( , , [ := = ( ) , 0 0 1 = ( ) , 0 1 2 = ( ) , 1 0 x = ( ) , 1 1 x2 ]) B[0,1]; 2 2. 向量定义 vector 格式: vector(n, [x1, x2, ..., xn]) vector(n, f) 说明: 参数中至少有向量维数 n 或初值表[x1, x2, ..., xn]. 95 vector(n)定义一个没有确定值的向量. vector(n, f)=vector(n,[f(1),f(2),...,f(n)]). 向量是下标从一开始的一维数组, 即vector(n)=array(1..n). 取向量 A 元素用 A[1], A[2]等等. vector 属于 linalg 软件包. 例with(linalg): #包含线性代数软件包 Warning, the protected names norm and trace have been redefined and unprotected vector([3,4,5,6,7]); [3, 4, 5, 6, 7] V:=vector(4, [1,3,5,7]); V:= [1, 3, 5, 7]; W:=vector(4,(i)->i^2); W:= [1, 4, 9, 16]; A:=array(1..4, [1,3,5,7]); A:= [1, 3, 5, 7]; type(V, array); #向量是一维数组 true type(A, vector); #下标从 1 开始的数组是向量 true 3. 矩阵定义 matrix 格式: matrix(m, n, [[a11,a12,...,a1n],[a21,a22,...,a2n]am1,am2,...,amn]]) matrix(m, n, f) 说明: 参数中至少有向量维数 m、n 或初值表[[a11,a12,...,a1n],[a21,a22,...,a2n]am1,am2,...,amn]] matrix(m, n)定义一个没有确定值的矩阵. 矩阵是下标从一开始的二维数组, 即matrix(m, n)=array(1..m, 1..n). 取矩阵 A 元素用 A[1, 1], A[1, 2]等等. matrix 属于 linalg 软件包. 例with(linalg): #包含线性代数软件包 Warning, the protected names norm and trace have been redefined and unprotected A:=matrix(3, 3); := A ( ) array , , .. 1 3 .. 1 3 [ ] B:=matrix([[1,2,3], [4,5,6], [7,8,9]]); := B ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 2 3 4 5 6 7 8 9 C:=matrix(3, 3, (i,j)->i^2+j^2); := C ? ? ? ? ? ? ? ? ? ? ? ? ? ? 2 5 10 5 8 13 10 13 18 type(C, matrix); #C 是矩阵类型 96 true type(C, array); #C 是二维数组 true 4. 矩阵的基本运算 matadd(A,B)或A+B —— 矩阵加法 multiply(A,B)或&*(A,B)或A&*B —— 矩阵乘法 scalarmul(A, λ)或λ*A —— 数乘矩阵 transpose(A) —— 转置矩阵 A inverse(A) —— 求A的逆矩阵 例with(linalg): A, B:=matrix([[1,2,3], [4,5,6], [7,8,9]]), matrix(3, 3, (i,j)->i^2+j^2); := , A B , ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 2 3 4 5 6 7 8 9 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 2 5 10 5 8 13 10 13 18 evalm(A+B); #矩阵加法 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 3 7 13 9 13 19 17 21 27 evalm(&*(A,B)); #矩阵乘法 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 42 60 90 93 138 213 144 216 336 transpose(A); #矩阵转置 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 4 7 2 5 8 3 6 9 inverse(matrix(2,2,[[a,b], [c,d]])); #矩阵求逆 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? d ? + a d b c b ? + a d b c c ? + a d b c ? a ? + a d b c 4. 矩阵的初等变换 addrow(A, r1, r2, m) —— 将矩阵 A 的r1 行乘上 m, 再加到 A 的r2 行上 mulrow(A, r, m) —— 将矩阵 A 的r行乘上 m swaprow(A, r1, r2) —— 将矩阵 A 的r1 与r2 行互换 注: 只要将函数名的"row"替换成"col"就可以进行相应的列变换. 5. 解线性联立方程组 AX = B 格式: linsolve(A, B) 求解解线性联立方程组 AX = B, 返回解向量, 属于 linalg 软件包. 例97 A:=matrix(3,3,[[1,0,2],[0,4,-3],[3,6,0]];B:=vector(3,[9,1,0]); 6. 矩阵的行列式、特征值和特征向量 det(A) —— 计算 A 的行列式值 charpoly(A) —— 求矩阵 A 的特征方程 eigenvalues(A) —— 求矩阵 A 的特征值 eigenvectors(A) —— 求矩阵 A 的特征向量 例A:=matrix([1,2,3],[1,2,3],[2,5,6]); 对eigenvectors 矩阵特征向量结果的注释. §1.5 微积分 1. 极限 limit 格式: limit(f(x), x=a, dir) —— 求)(lim x f a x→ , dir 指定逼近方向 说明: dir 为逼近方向, 其值可以是 left(左极限)、 right(右极限)、 real(双方向极限)与complex(全 方向极限). 缺省则以双方向逼近. 逼近的点可以是一个常数、变量或者是 infinity(∞)、-infinity(∞). complex 用在多元函数 Maple 允许极限为 infinity(∞) 例limit(sin(x)/x, x=0); 1 limit(ln(abs(x))/x, x=0); #极限不存在 undefined limit(ln(abs(x))/x, x=0,right); #右极限为无穷大 -∞ limit((x+1)/x, x=infinity); 1 limit((sin(x)-sin(a))/(x-a), x=a); cos(a) a:=(n)->1/sqrt(5)*(((1+sqrt(5))/2)^n-((1-sqrt(5))/2)^n); #费波那奇数列通项公式 := a → n ? ? ? ? ? ? ? ? ? + 1 2 1 2 5 n ? ? ? ? ? ? ? ? ? 1 2 1 2 5 n 5 limit(a(n)/a(n+1), n=infinity); 2 1 + 1 5 98 evalf(%); #黄金分割点 .6180339886 2. 微分 diff、D 格式: diff(f(x),x) —— 求()nndfxdx , diff(f(x),x$n) —— 求()nndfxdx D(f) —— 计算函数 f 的一阶微分 (D@@n)(f) —— 计算函数 f 的n阶微分和 说明: diff(f(x),x)是对表达式求微分, 而D(f)中的 f 是函数名, 不是 f(x). diff(f(x1,x2,...),x1,x2,...)和D[1,2,...](f)用来求多元函数微分 例diff(sin(x), x); cos(x) diff(sin(x), x$2); -sin(x) f:=D(sin); f:=cos g:=(D@@2)(sin); g:=-sin f(Pi); -1 diff('f'(x,y), x, y); ? ? ?2 y x ( ) f , x y diff(sin(x*y), x, y); ? + ( ) sin x y x y ( ) cos x y s:=(x,y)->sin(x*y); := s → ( ) , x y ( ) sin x y D[1,2](s); → ( ) , x y ? + ( ) sin x y y x ( ) cos x y 3. 积分 格式: int(f(x), x) —— 求不定积分 ∫ dx x f ) ( int(f(x), x=a..b) —— 求定积分 ∫ b a dx x f ) ( Doubleint(f(x,y), x, y) —— 求二重不定积分 ∫∫ dxdy y x f ) , ( Doubleint(f(x,y), x=a..b, y=c..d) —— 求二重定积分 ∫ ∫ d c b a dy dx y x f ) , ( 99 说明: int(f(x), x)计算不定积分时 Maple 省去了积分常数 C. Doubleint()只求出积分形式, 并没有 作真正的计算, 可用 value()求出值. int 与Doubleint 属于 student 软件包. 例with(student): int(1/(x^3+1), x); ? + 1 3 ( ) ln + x 1 1 6 ( ) ln ? + x2 x 1 1 3 3 ? ? ? ? ? ? ? ? arctan 1 3 ( ) ? 2 x 1 3 int(1/(x^3+1), x=0..1); + 1 3 ( ) ln 2 1 9 3 π Doubleint(f(x,y), x, y); d ? ? ? ? d ? ? ? ? ( ) f , x y x y Doubleint(x*y, x, y); d ? ? ? ? d ? ? ? ?x y x y value(%); 1 4 x2 y2 5. 微分方程求解 格式: dsolve({eqns, init}, f(x)) —— 求解微分方程 eqns, init 是初始条件, f(x)指定要求解的函数 pdsolve(eqns, f(x)) —— 求解偏微分方程 eqns, f(x)指定要求解的函数 例eqn:=diff(y(x),x$2)+2*diff(y(x),x)+1=0; := eqn = + + ? ? ? ? ? ? ? ? ? ? ? ?2 x2 ( ) y x 2 ? ? ? ? ? ? ? ? ? ? x ( ) y x 1 0 dsolve(eqn); = ( ) y x ? ? + 1 2 e ( ) ?2 x _C1 1 2 x _C2 dsolve({eqn, y(0)=0, (D@@2)(y)(0)=1}); #解含初值条件的微分方程 = ( ) y x ? ? 1 4 e ( ) ?2 x 1 2 x 1 4 PDE:=x*diff(f(x,y),y)-diff(f(x,y),x)=f(x,y)^2*g(x)/h(y); := PDE = ? x ? ? ? ? ? ? ? ? ? ? y ( ) f , x y ? ? ? ? ? ? ? ? ? ? x ( ) f , x y ( ) f , x y 2 ( ) g x ( ) h y pdsolve(PDE); #解偏微分方程 = ( ) f , x y 1 + d ? ? ? ? ? ? ? ? ? x ( ) g _a ? ? ? ? ? ? ? ? h ? + + 1 2 _a2 y 1 2 x2 _a ? ? ? ? ? ? ? ? _F1 + y 1 2 x2 100 §1.6 作图 1. 二维曲线作图 plot 格式: plot(f(x), x=a..b) —— 作f(x)的图形, x 在[a, b]内plot([f1(x), f2(x)x=a..b) 在同一坐标系中作多个函数图形 plot([x(t), y(t), t=a..b]) —— t 从a到b作二维参数图 plot([[x1,y1], [x2,y2]style=point) —— 数据点绘图 plot(f(x), x=a..b, options) —— 加入选项来更改绘图指令的默认值 说明: options 的形式如 option=value, 取值如下 —— axes 指定坐标轴形状: FRAME(框), BOXED(盒子), NORMAL(普通), NONE(无) —— color 指定作图颜色, aquamarine(碧绿), black(黑), blue(蓝), navy(藏青), coral(珊瑚), cyan(墨绿), brown(棕色), gold(金), green(绿), gray(灰白), grey(灰), khaki(土黄), magenta(洋红), maroon(栗色), orange(橘黄), pink(粉红), plum(梅红), red(红), sienna(赫色), tan(黄褐), turquoise(蓝玉), violet(紫), wheat(灰黄), white(白), yellow(黄) —— coords 坐标系选 择, cartesian(笛卡 尔坐 标系, 缺省), elliptic(椭圆坐标系 ), logarithmic(对数坐标系), polar(极坐标) 详细的说明请在 Maple 中键入?plot[options]来查询 options 不仅可以用在 plot 函数中, 几乎所有的 Maple 作图函数都可以使用 例plot(sin(x),x=0..2*Pi); plot([sin(x)+cos(x)^2, (sin(x)+cos(x)^2)^2, (sin(x)+cos(x)^2)^3, (sin(x)+cos(x)^2)^4, (sin(x)+cos(x)^2)^5], x=0..16); data:=[1,1],[2,4],[3,9],[4,16],[5,25],[6,36],[7,49],[8,64],[9,81],[10,100]: plot([data],style=point,color=blue); 101 plot([seq(sin(5*t)+n*cos(t),n=-5..5)],t=0..Pi,coords=polar); 2. 三维作图 plot3d 格式: plot3d(f(x,y), x=a..b, y=c..d) —— 在[a, b]*[c, d]内作 f(x,y)的图形 plot3d(f(x,y), x=a..b, y=c..d, options) —— 加入选项来更改绘图指令的默认值 说明: 用法可参考 plot options 取值如下 —— ambientlight=[r, g, b] 设置光线强度, 数值范围是 0..1 —— axes 同plot —— color 同plot —— coords 指定坐标系, spherical(球坐标), cylindrical(圆柱坐标) 详细的说明请在 Maple 系统中键入?plot3d[options]查询 例plot3d(sin(x)*cos(y),x=-Pi..Pi,y=-Pi..Pi); plot3d(sin(t)*sin(p^2),t=0..Pi,p=0..Pi,coords=spherical); #球坐标作图 102 plot3d(sin(z/2),t=0..2*Pi/2,z=-Pi..Pi,coords=cylindrical); #圆柱坐标作图 3. 显示图形 display 格式: display(g1, g2,显示图形 g1, g2, ... 在同一坐标下 说明: 如果多个图形无法一次画出, 但又想将它们显示在同一坐标下, 可以使用 display display 既可以显示二维图, 也可以显示三维图, 但不能将二维图与三维图同时显示 display 属于 plots 软件包 例with(plots): g1:=plot3d(2*exp(-sqrt(x^2+y^2)), x=-6..6,y=-6..6): g2:=plot3d(sin(sqrt(x^2+y^2)), x=-6..6,y=-6..6): display(g1, g2); 4. Maple 动画 103 格式: animatecurve(f(x), x=a..b) —— 模拟手绘函数 f(x)的过程 animate(f(x,t), x=a..b, t=t1..t2) —— 绘制二维平面动画, t 代表时间 animate3d(f(x,y,t), x=a..b, y=c..d, t=t1..t2) —— 绘制三维动画, t 代表时间 说明: 同属于 plots 软件包. 例with(plots): animatecurve([sin(t),cos(t),t=-Pi..Pi]); animate(sin(x*t), x=0..2*Pi, t=0..8, coords=polar, frames=140); animate3d(cos(x)*sin(y)*sin(t), x=-Pi..Pi, y=-Pi..Pi, t=0..2*Pi, frames=36); (图略) 运行动画函数后 Maple 并没有马上显示动画, 如果想观看动画, 须用鼠标选择图形, 此时在工具 栏上出现播放动画控制按扭 , 点击即可播放. 5. 图形的保存 Maple 可以把图像转换成各种图形文件格式保存起来. 方法: 用鼠标右键单击图形, 然后选择 Export As, 选择要转换的格式, 再在对话框中选择好保存路径 并键入文件名, 点击确定即可. 如果图形是动画, 并且保存为 GIF 格式, 就可以将 Maple 动画用于网页制作及其它应用了. §1.7 函数与过程 Maple 的函数与过程其实并没有严格的区别, 一般来说, 函数用"->"箭头操作符来定义, 而过程 用关键字 proc 来定义. Maple 的函数与过程的调用方法是完全一样的. a)函数 1. 函数定义 格式: f := (x, y,expr —— 定义函数 f(x, y, ... ) 例f:=(x)->x^2+2*x-1; := f → x + ? x2 2 x 1 f(sin(x)); + ? ( ) sin x 2 2 ( ) sin x 1 plot(f(x),x=-3..3); g:=(x,y)->sin(x)+cos(y)^2; 104 := g → ( ) , x y + ( ) sin x ( ) cos y 2 plot3d(g(x,y),x=-Pi..Pi,y=-Pi..Pi); b)过程 1. 过程定义 格式: p := proc( P ) —— 过程定义语句, proc 是关键字, P 为形式参数列表 local v1, v2,定义局部变量 global g1, g2,定义全局变量 option o1, o2,关于过程的选项, 用于过程记忆表, 操作符, 内部函数组成, 执 行跟踪, 版本说明等 proc body —— 过程体 end proc; —— 过程定义结束语句 说明: 由于过程定义较为复杂, 在此将各部分分开解释 (1) 形式参数 格式: p := proc(parem_1::type_1, ..., parem_k::type_1) —— parem 变量, type 类型 proc body end proc; 说明: 调用过程 p(e_1, ..., e_2)的时候, 实际参数表的值从左到右依次传递给形式参数表, 即每 一个 e_i 的值传给 parem_i, 然后检测每个参数的类型. 只要其中之一不匹配, 则返回出错信息. 常用的参数类型 type 有: numeric(数值), float(浮点数), integer(整数), posint(正整数), positive(正数), relacons(实数), complex(复数), even(偶数), odd(奇数), equation(方程), boolean(布尔)等, 更详细的请运 行"?type"查看帮助. 参数在过程中不能再赋值. 实际参数的数目可以与形式参数的数目不同, 可以多于但不得少于. 不需要宣告参数类型时, 用形式 proc(x)与proc(x::anything)相同. 使用变量名 args, 允许用户访问过程内部的实际参数, args[1]取得第一个参数值, args[2]取得第二 参数值, 以此类推. 也就是说, 允许用户不定义形式参数, 但仍然可以访问实际参数. 使用变量 nargs, 可以得到实际参数个数. 因此, 使用变量 args 和nargs, 程序变换变得更为灵活 例f:= proc(x, y, z) x+y+z; end proc; := f proc ( ) end proc , , x y z + + x y z 105 f(a, b, c); f(1, 2, 3, 4); #实参个数允许多于形参个数, 但此时只有前 3 个有效 6 MAX := proc(x::numeric, y::numeric) if x>y then x; else y; end if; end proc; MAX(Pi, 3); #出错, 因为 Pi 不是 numeric 类型 Error, invalid input: MAX expects its 1st argument, x, to be of type numeric, but received Pi MAX(evalf(Pi), 3); #将Pi 转化为浮点数后再做比较 3.141592654 例3#使用变量 args 和nargs, 不提供形式参数而直接使用实际参数 #求最大值的函数, 参数个数可以是有限个 MAX := proc() local i, m; if nargs=0 then return FALL end if; m:=args[1]; for i from 2 to nargs do if args[i]>m then m:=args[i] end if; end do; return m; end proc: MAX(evalf(Pi), 22/7, 104348/33215); MAX(evalf(Pi), 22/7, 104348/33215, 312689/99532, 5419351/1725033); (2) 全局变量和局部变量 变量在过程中要么是局部变量, 要么是全局变量. Maple 对于在不同的过程里有相同的变量名的 局部变量, 视为不同的变量. 这样, 在一个过程里的局部变量值的改变, 不影响其它同名的变量值. 格式: p := proc( ) local v1, v2,定义局部变量 global g1, g2,定义全局变量 proc body end proc; Maple 允许用户不申明变量而直接使用, 那么在过程中未经申明的变量是全局变量还是局部变 量? Maple 自动将下列两种情形视为局部变量, 否则视为全局变量: —— 变量出现在赋值号的左边 —— 变量作为下标变量出现在 for 循环中, 或seq, add, mul 等函数中 (3) 过程选项 格式: p := proc( ) 106 option o1, o2,选项, 用于过程记忆表, 操作符, 内部函数组成, 执行跟踪, 版 本说明等 proc body end proc; 说明: 选项有: remember(记忆表), builtin(内部函数组成), system(系统管理), operator(操作符), arrow(箭头符), trace(执行跟踪), Copyright(版本说明). 在此只解释记忆表 remember, 如想了解其它的请自行查阅帮助文档. 使用过程选项 remember, Maple 就在内存中建立一个表(表是类似数组的一种数据结构), 用来存 储过程的返回值, 每遇到一个新的实参, 就把这个参数作为表的索引(下标), 把过程的返回值作为表 的值. Maple 每执行一次这个过程, 就检查程序是非否存在以前程序执行的结果. 如果程序已经保存 这些计算结果, Maple就直接取回结果而不是再执行同一参数的过程. 这对于执行递归过程特别有用, 使之效率倍增. 使用 op(4, eval(p))可以观察过程 p 的记忆表的结构, 甚至可对记忆表进行各种操作, 操作方法请 查阅相关帮助文档. 例 #定义费泊那契数列, 使用递归定义 fib := proc(n::nonnegint) #参数 n, 限定非负整数 option remember; #使用 remember 选项 if n<2 then n else fib(n-1) + fib(n-2) end if; end: fib(1); 1 fib(4); 3 fib(6); 8 op(4, eval(fib)); #显示记忆表 table([0 = 0, 1 = 1, 2 = 1, 3 = 2, 4 = 3, 5 = 5, 6 = 8]) fib(7); 13 fib(8); 21 op(4, eval(fib)); #显示记忆表 table([0 = 0, 1 = 1, 2 = 1, 3 = 2, 4 = 3, 5 = 5, 6 = 8, 7 = 13, 8 = 21]) fib(4); #当记忆表中有参数为 4 的值时, 过程内的语句不被执行, 而是直接从记忆表中取出结果 3 (4) 过程返回值 格式: RETURN(v_1, v_2, ..., v_n) —— 返回一个或多个值 v_1, v_2, ..., v_n return v_1, v_2, ..., v_n —— 返回一个或多个值 v_1, v_2, ..., v_n ERROR(v_1, v_2, ..., v_n) —— 出错信息 WARNING(msg)—— 警告信息 说明: RETURN 的普遍形式是返回过程执行的最后一条语句, RETURN 与return 无多大差别. 用RETURN('procname'(args))可以实现返回过程调用原样. 返回多个值时是返回一个序列, 如果执行 v:=p()而p()有多个返回值, 那么 v[1], v[2]分别 107 取得第 1、2 个返回值. ERROR 和WARNING 用于出错处理. §1.8 流程控制 1. 条件判断语句 if 格式: if cond_1 then statement_1; elif cond_2 then statement_2; ... ... else else_statement; end if; 说明: cond_i 是布尔表达式, statement_i 是语句组. 当cond_1 为真(true)时, 执行语句 statement_1, 否则当 cond_2 为真时, 执行语句 statement_2, ..., 当所有的布尔表达式都为假(false)时, 执行语句 else_statement. elif 和else 都是可选项, 如果不须要的话可以省去. 一种简化的 if 形式: `if`(cond, x, y) —— 如果 cond 为真(true), 返回 x, 否则返回 y. 注意这里的 if 必须用反单引号(在 键盘左上角 1 旁边的键)括起来. 例ABS := (x) -> if (x>0) then x; else -x; end if; #绝对值函数 ABS(3); 3 ABS(-3); 3 ABS(a); #当参数为符号时出错 Error, (in ABS) cannot evaluate boolean: -a < 0 ABS := (x) -> if not type(x, numeric) then #改进的绝对值函数 RETURN('procname'(args))(x>0) else if (x>0) then x; else -x; end if; end if; ABS := proc (x) options operator, arrow; if not type(x,numeric) then RETURN(('procname')(args))(0 < x) else if 0 < x then x else -x end if end if end proc ABS(-4); 4 ABS(a); ABS(a) 2. 循环语句 格式: for var from v0 by step to v1 while cond do 108 statement end do; 说明: 变量 var 从v0 到v1, 步长为 step, 条件为 cond, 以循环的方式执行 statement 除了 do statement end do 外其它语句都是可选项, for 缺省为虚拟变量, from 和by 缺省为 1, to 缺 省为 infinity(∞), while 缺省为 true. 例for i from 10 by 5 to 30 do print(i) end do; 10 15 20 25 30 i:=5; i:=5 while i<20 do i:=i+5; end do; i:=10 i:=15 i:=20 #一个求<=n 的所有质数的例子 prim := proc(n) local p, i; p[1]:=2; for i do p[i+1]:=nextprime(p[i]); if nextprime(p[i+1])>n then break; end if; end do; return op(p); end proc: p := prim(13); p[1]; 2 p[6]; 13 当遇到了 break 命令, 退出当前循环. 第三节 LINDO 与LINGO §1.1 LINDO 与LINGO 简介 LINDO 是Linear INteractive and Discrete Optimizer 字首的缩写形式, 是由 Linus Schrage 于1986 年开发的优化计算软件包, 可以用来求解线性规划 (LP——Linear Programming), 整数规划 (IP——Integer Programming) 和二次规划 (QP——Quadratic Programming) 问题. LINDO 易于规划问 109 题的输入、 求解和分析, 程序执行速度很快. LINDO 学生版最多可求解多达 200 个变量和 100 个约束 的规划问题. LINGO 可用于求解线性规划和整数规划问题. LINGO 与LINDO 不同的是, LINGO 包含了内置的建模语言, 允许以简练、直观的方式描述所 需求解的问题, 模型中所需的数据可以以一定格式保存在列表(List)和表格(Table)中, 也可以保存在 独立的文件中. §1.2 用LINDO 求解线性规划(LP)问题 初试 LINDO LINDO 的求解机制: LINDO 的求解过程采用单纯形法, 一般是首先寻求一个可行解, 在有可 行解情况下再寻求最优解. 用LINDO 求解一个 LP 问题会得到如下的几种结果: 不可行(No feasible solution) 或 可行(Feasible) 可行时又可分为: 有最优解(Optimal Solution)和解无界(Unbounded Solution)两种情况. 由于在实 际问题中, 不太可能出现最大利润无上限的情形, 所以使用者应检查是否少了一个约束或有其它印 刷错误. 让我们来解如下 LP 问题: max z=2x+3y s.t. 4x+3y≤10 3x+5y≤12 x,y≥0 由于 LINDO 中已假设所有的变量都是非负的, 所以非负约束可不必再输入到计算机中; LINDO 也不区分变量中的大小写字符(实际上任何小写字符将被转换为大写字符); 约束条件中的"<=" 及">="可用"<" 及">"代替. 上面问题用键盘输入如下: MAX 2X + 3Y ST 4X + 3Y < 10 3x + 5Y < 12 END LINDO 中一般称上面这种问题实例(INSTANCE)为模型(MODEL). 以后涉及该模型时, 目标函 数为第一行, 两个约束条件分别为第二、三行. 点击 求解模型, 结果如下: LP OPTIMUM FOUND AT STEP 2 OBJECTIVE FUNCTION VALUE 1) 7.454545 VARIABLE VALUE REDUCED COST X 1.272727 0.000000 Y 1.636364 0.000000 ROW SLACK OR SURPLUS DUAL PRICES 2) 0.000000 0.090909 3) 0.000000 0.545455 NO. ITERATIONS= 2 计算结果表明: "LP OPTIMUM FOUND AT STEP 2"表示单纯形法在两次迭代(旋转)后得到最优解. 110 "OBJECTIVE FUNCTION VALUE 1) 7.4545450 "表示最优目标值为 7.4545450. "VALUE"给出最优解中各变量(VARIABLE)的值: X =1.272727, Y =1.636364. "REDUCED COST" 给出最优单纯形表中第 0 行中变量的系数 ( max 型问题). 其中基变量的 reduced cost 值应为 0, 对于非基变量, 相应的 reduced cost 值表示当该非基变量增加一个单位时目标 函数减少的量. 本例中此值均为 0. "SLACK OR SURPLUS" 给出松驰变量的值: 第2、3 行松驰变量均为 0, 说明对于最优解来讲, 两个约束(第2、3 行)均取等号. "DUAL PRICES" 给出对偶价格的值: 第2、3 行对偶价格分别为 .090909, .545455. "NO. ITERATIONS= 2" 表示用单纯形法进行了两次迭代(旋转). 一个问题解答之后, LINDO 会询问是否需要做灵敏性分析(DO RANGE (SENSITIVITY) ANALYSIS? ) 如果你不需要, 你应回答"N"(NO). 存储模型 下面即是一个具体应用的例子: (可参照上述使用步骤) 首先输入问题: MAX 60 DESKS + 30 TABLES + 20 CHAIRS SUBJECT TO 2) 8 DESKS + 6 TABLES + CHAIRS <= 48 3) 4 DESKS + 2 TABLES + 1.5 CHAIRS <= 20 4) 2 DESKS + 1.5 TABLES + 0.5 CHAIRS <= 8 5) TABLES <= 5 END 问题求解, 要求做敏感性分析 LP OPTIMUM FOUND AT STEP 2 OBJECTIVE FUNCTION VALUE 1) 280.0000 VARIABLE VALUE REDUCED COST DESKS 2.000000 0.000000 TABLES 0.000000 5.000000 CHAIRS 8.000000 0.000000 ROW SLACK OR SURPLUS DUAL PRICES 2) 24.000000 0.000000 3) 0.000000 10.000000 4) 0.000000 10.000000 5) 5.000000 0.000000 NO. ITERATIONS= 2 RANGES IN WHICH THE BASIS IS UNCHANGED: OBJ COEFFICIENT RANGES VARIABLE CURRENT ALLOWABLE ALLOWABLE COEF INCREASE DECREASE DESKS 60.000000 20.000000 4.000000 TABLES 30.000000 5.000000 INFINITY CHAIRS 20.000000 2.500000 5.000000 RIGHTHAND SIDE RANGES ROW CURRENT ALLOWABLE ALLOWABLE 111 RHS INCREASE DECREASE 2 48.000000 INFINITY 24.000000 3 20.000000 4.000000 4.000000 4 8.000000 2.000000 1.333333 5 5.000000 INFINITY 5.000000 仍以上面的问题 DAKOTA 为例, 下面给出其结果的一般注释: "LP OPTIMUM FOUND AT STEP2"表示 LINDO 在(用单纯形法)两次迭代或旋转后得到最优解. "OBJECTIVE FUNCTION VALUE 280.000000"表示最优目标值为 280. "VALUE"给出最优解中各变量的值. 例Dakota 问题中需造 2 个(书桌)desks, 0 个(桌子)tables, 和8个(椅子) chairs. "SLACK OR SURPLUS"给出松驰变量的值. 上例中 s1= 第2行松驰变量 =24 s2= 第3行松驰变量 =0 s3= 第4行松驰变量 =0 s4= 第5行松驰变量 =5 "REDUCED COST"出最优单纯形表中第 0 行中变量的系数(max 型问题). 其中基变量的 reduced cost 值应为 0, 对于非基变量, 相应的 reduced cost 值表示当增加一个单位时目标函数减少的量. 敏感性分析 使用 LINDO 时, 结果输出中会提供敏感性分析. 这一信息一般包含于两个标题之下, 其一是 REDUCED COSTS, 另一个是 DUAL PRICES. 它们分别表示了当变量或约束条件有微小变动时, 目 标函数的变化率. 在输出结果中对应于每个变量都有一个 REDUCED COST, 若其数值为 x, 表示对应的变量为零 时, 若增加 1 个单位, 目标函数将减少 x 个单位. 输出结果中对应于每一个约束也都有一个 DUAL PRICE. 若其数值为 x, 表示对应约束中不等式 右端项若减少 1 个单位, 目标函数将增加 x 个单位. 如果 REDUCED COST 或DUAL PRICE 的值为 0, 表示微小扰动不影响目标函数. 有时, 通过分析 DUAL PRICE, 也可对产生不可行问题的原因有所了解. 注意事项: 1) LINDO 中已假定所有变量非负. 变量名不能超过 8 个字符. 2) 如要输入 <= 或>= 型约束, 相应以< 或 >代替即可. 3) LINDO 不允许变量出现在一个约束条件的右端. 4) 目标函数及各约束条件之间一定要有空格分开. 5) 一般 LINDO 中不能接受括号( )和逗号"," , 例:400(X1+X2)需写为 400X1+400X2; 10,000 需 写为 10000. 6) 14)LINDO 将目标函数所在行作为第一行, 从第二行起为约束条件. 行号自动产生, 也可以 人为定义行号或行名. 行名和变量名一样, 不能超过 8 个字符. 7) 15)数值均衡化及其它考虑. LINDO 不能将 LP 中的矩阵进行数值均衡化. 为了避免数值问 题, 使用者应自己对矩阵的行列进行均衡化. 一个原则是, 系数矩阵中非零元的绝对值不能大于 100,000 或者小于 0.0001. 如果 LINDO 觉得矩阵元素之间很不均衡, 将会给出警告. 8) 量纲分析与一般错误的避免. LINDO 部分常命令中文注释如下: MAX/MIN 命令: 用于输入一个包含目标函数,约束条件在内的 LP 模型. 输入程序如下: 输入"MAX" (或"MIN"), 继之以自然格式的目标函数作为第一行;再输入 "SUBJECT TO"(可简写为"ST"), 后面跟约束条件行. 最后, 输入"END"回到命令状态模式. 以后只需 112 给出"GO"命令即可开始优化求解过程. 其中, 变量名可以由 1—8 个字母或数字型的字符构成, 且第一个字符必须是字母. 变量系数不 能是指数型, 例如: .258E+29 形式的系数是不允许的. 关键词 ("MAX","ST","END"...)及各行之间必 须用一个或多个空格分隔开. 空格可以出现在一行之中, 但不能出现在变量名中. 一个回车符等价 于一个空格. 下面是同一问题的两种合法的输入方式: 1) MIN 2X+3Y SUBJECT TO -5X-2Z<=10 +10X - Y >5 END 2) MIN 2X + 3 Y ST -5X-2Z<10 10X -Y>+5 END 其它一些有用的符号有: 算术运算符 逻辑运算符 关系运算符 顺序运算符 AND. LOG( ) EXP( ) .OR. ABS( ) .NOT. GIN 命令: GIN 命令可将问题模型中的变量标为 0/1 型. 第一种格式为"GIN n" , 其中 n 是整型变量的个数. LINDO 要求整型变量应放在问题模型的最前面. 第二种格式为 "GIN var-id", 其中"var-id"是变量名. §1.3 用LINDO 求解整数规划和二次规划问题 3.1 整数规划(IP) LINDO 可用于求解单纯的或混合型的整数规划(IP)问题. 但目前尚无相应完善的敏感性分析理 论. IP 问题的输入与 LP 问题类似, 但在 END 标志后需定义整型变量. 0/1 型的变量可由 INTEGER(可简写为 INT)命令来标识: INTEGER vname 或INTEGER n 前者只将变量 vname 标识为 0/1 型, 后者将当前模型中前 n 个变量标识为 0/1 型. 模型中变量顺 序由输入决定, 该顺序可由输出结果中查证. 具体演示如下: MAX 4TOM +3DICK +2HARRY ST 2.5TOM +3.1HARRRY < 5 .2TOM +.7DICK +.4HARRY < 1 END :INT TOM :INT DICK :INT HARRY 求解结果: LP OPTIMUM FOUND AT STEP 3 OBJECTIVE VALUE = 7.71428585 113 NEW INTEGER SOLUTION OF 7.00000000 AT BRANCH 0 PIVOT 3 RE-INSTALLING BEST SOLUTION... OBJECTIVE FUNCTION VALUE 1) 7.000000 VARIABLE VALUE REDUCED COST TOM 1.000000 -4.000000 DICK 1.000000 -3.000000 HARRY 0.000000 -2.000000 HARRRY 0.000000 0.000000 ROW SLACK OR SURPLUS DUAL PRICES 2) 2.500000 0.000000 3) 0.100000 0.000000 NO. ITERATIONS= 3 BRANCHES= 0 DETERM.= 1.000E 0 注: 整数规划中变量不一定限于 0/1 型, 因此 LINDO 中有一个比 INT 弱些的命令 GIN, 可将变 量仅限为整数型, 而其使用方式及格式与 INT 命令相似. 尽管 LINDO 对整数规划问题是很有威力的, 要想有效地使用还是需要一定的技术. 这是因为, 人们很容易将一个本质上很简单的问题列成一个不好的数学表达式, 而一个坏的表达式有可能会导 致一个冗长的计算. 当然这时 LINDO 会主动砍去一些过程, 以缩短计算时间. 3.2 二次规划(Quadratic Programming) Lindo 可用于求解二次规划问题. 例如: 对于 min z=3x2 +y2 -xy+0.4y s.t. 1.2x+0.9y>1.1 x+y=1 y<0.7 该问题通过在实际约束前增加有关变量的一阶条件, 从而转化二次型为线性(互补)型(如下). 这 需要我们为每一个实际约束增加一个对偶变量, 并要使用 QCP 命令. 第一行(目标函数)只用于给出 相应变量的顺序. 对上面的问题, 我们将用 RT, ONE 和UL 作为对偶变量. 问题输入如下: MIN X+Y+RT+ONE+UL ST 6X - Y - 1.2RT + ONE > 0. - X + 2Y - .9RT + ONE + UL > -.4 1.2X + .9Y > 1.1 X + Y = 1 Y < .7 END 求解: LP OPTIMUM FOUND AT STEP 1 OBJECTIVE FUNCTION VALUE 1) 1.000000 VARIABLE VALUE REDUCED COST X 0.666667 0.000000 Y 0.333333 0.000000 RT 0.000000 1.000000 114 ONE 0.000000 1.000000 UL 0.000000 1.000000 ROW SLACK OR SURPLUS DUAL PRICES 2) 3.666667 0.000000 3) 0.400000 0.000000 4) 0.000000 0.000000 5) 0.000000 -1.000000 6) 0.366667 0.000000 NO. ITERATIONS= 1 §1.4 LINGO 功能简介 LINGO 软件包的功能之一在于提供了一些复杂的矩阵生成器的例子, 更进一步说, LINGO 实际 上提供了建立最优化模型的一种语言, 有了它, 使用者只用键入一行文字也可以建立起成千条约束 或目标函数项. 这就使输入较大规模问题的过程得到了简化. 一般地, LINGO 的内部函数有: 普通函数 @ABS( X) @BIN( VAR) @SIN( X) @COS( X) @EXP( X) @FPA( I, N) @FPL( I, N) @FREE( VAR) @GIN( VAR) @IMPORT( WK1 FILE, RANGE NAME) @IN( SET NAME, SET ELEMENT) @LGM( N) @LOG( X) @PBN( P, N, X) @PEB( ARL, S) @PEL( ARL, S) @PHG( POP, GREEN, SAM, X) @PFS( ARL, S, C) @PPL( ARL, S) @PPS( ARL, S) @PSN( Z) @PSL( Z) @SIZE( SET) @SMAX( X, Y) @SMIN( X, Y) @TAN( X) @USER( X1, ..., XN) 115 @WARN( 'TEXT', CONDITION) @WRAP( INDEX, SET SIZE) 集合函数 @MAX( SET | COND: EXP) @MIN( SET | COND: EXP) @SUM( SET | COND: EXP) @FOR( SET | COND: EXP) 例: LINGO 中问题的一个典型模型如下: (SAILCO) 4 1 1 1 1 1 1 MIN 400 450 20 S.T. 40, 1..4 , 1, 1..4 10 i i i i i i i i i i i RP OP INV RP i INV INV RP OP DEM TIME i INV RP OP DEM = ? + + < = = + + ? ∑ 其中, DEM={40,60,75,25}, TIME={1,2,3,4}. 对应的 LINGO 模型如下: MODEL: SETS: QUARTERS/Q1,Q2,Q3,Q4/:TIME,DEM,RP,OP,INV; ENDSETS MIN=@SUM(QUARTERS:400*RP+450*OP+20*INV); @FOR(QUARTERS(I):RP(I)<40); @FOR(QUARTERS(I)|TIME(I)#GT#1: INV(I)=INV(I-1)+RP(I)+OP(I)-DEM(I);); INV(1)=10+RP(1)+OP(1)-DEM(1); DATA: DEM=40,60,75,25; TIME=1,2,3,4; ENDDATA END 从上述 MODEL 中, 我们可看到整个模型分为三大部分: 1)SETSENDSETS 其中定义了模型中用到的各个集合, 包括变量, 数组, 变量的特征量等. 如上例中定义了四个 quarters:Q1,Q2,Q3,Q4. 其中每个 quarter 都有 TIME,DEM,RP,OP,INV 这样的特征量. 一旦这样的定 义建立起来, 实际上 quarter 的数量可以为 40, 400, 或许 4000, 它们仍都有 TIME, DEM, RP, OP, INV 这样的特征量. 这些量的具体数值可在 DATAENDDATA 这部分输入. 2)中间部分: 这部分实际上定义了目标函数, 约束条件等. 如上例中: 第4行定义了目标函数为 MIN=..., 其中@SUM(quarters: ... 400*RP...)表示对所有的 quarters 计算 400*RP+450*OP+20*INV 并求和. 注意, 如此定义的目标函数与 quarters 的数目是 4, 40, 400, 或4000 并无具体的关系. 第5行表示对每个 quarter, RP 不能超过 40 . 第6, 7 行也是对每个 quarter 中具体约束的定义. 116 注意: 与LINDO 不同的是变量可以放在约束条件的右端(同时数字也可放在约束条件的左端). 3)DATAENDDATA 行9-12 的作用在于输入必要的数据. 这部分要以 DATA: 开始, 以ENDDATA 结束. 求解结果如下: LP OPTIMUM FOUND AT STEP 7 OBJECTIVE VALUE = 78450.0000 VARIABLE VALUE REDUCED COST TIME( Q1) 1.000000 .0000000 TIME( Q2) 2.000000 .0000000 TIME( Q3) 3.000000 .0000000 TIME( Q4) 4.000000 .0000000 DEM( Q1) 40.00000 .0000000 DEM( Q2) 60.00000 .0000000 DEM( Q3) 75.00000 .0000000 DEM( Q4) 25.00000 .0000000 RP( Q1) 40.00000 .0000000 RP( Q2) 40.00000 .0000000 RP( Q3) 40.00000 .0000000 RP( Q4) 25.00000 .0000000 OP( Q1) .0000000 20.00000 OP( Q2) 10.00000 .0000000 OP( Q3) 35.00000 .0000000 OP( Q4) .0000000 50.00000 INV( Q1) 10.00000 .0000000 INV( Q2) .0000000 20.00000 INV( Q3) .0000000 70.00000 INV( Q4) .0000000 420.0000 ROW SLACK OR SURPLUS DUAL PRICE 1 78450.00 1.000000 2 .0000000 30.00000 3 .0000000 50.00000 4 .0000000 50.00000 5 15.00000 .0000000 6 .0000000 450.0000 7 .0000000 450.0000 8 .0000000 400.0000 9 .0000000 430.0000 LINGO 使用注意事项: 1) 所有的语句除 SETS, ENDSETS, DATA , ENDDATA, 和END 之外必须以一个分号 ";" 结尾. 2) 与LINDO 不同的是变量可以放在约束条件的右端(同时数字也可放在约束条件的左端). 3) 记住用一个分号";"将目标函数与约束条件分开. 4) LINGO可用于求解整数规划问题. 注意在MODEL中定义一个0-1型变量用 @BIN算子, 定 义一个非负整数用 @GIN 算子. 5) LINGO 解非线性规划时已假定各变量非负.