第1篇 建立数学模型
随着科学技术的迅速发展，数学模型这个词汇越来越多地出现在现代人的生产、工作和社会活动中。电气工程师必须建立所要控制的生产过程的数学模型，用这个模型对控制装置作出相应的设计和计算，才能实现有效的过程控制。气象工作者为了得到准确的天气预报，一刻也离不开根据气象站、气象卫星汇集的气压、雨量、风速等资料建立的数学模型。生理医学专家有了药物浓度在人体内随时间和空间的变化的数学模型，就可以分析药物的疗效，有效地指导临床用药。城市规划工作者需要建立一个包括人口、经济、交通、环境等大系统的数学模型，为领导层对城市发展规划的决策提供科学根据。厂长经理们要是能够根据产品的需求状况、生产条件和成本、贮存费用等信息，筹划出一个合理安排生产和销售的数学模型，一定可以获得更大的经济效益。就是在日常活动如访友、采购当中，人们也会谈论找一个数学模型，优化一下出行的路线。对于广大的科学技术人员和应用数学工作者来说，建立数学模型是沟通摆在面前的实际问题与他们掌握的数学工具之间联系的一座必不可少的桥梁。
本章作为全书的导言和数学模型的概述，主要讨论建立数学模型的意义、方法和步骤，给读者以建立数学模型的全面的、初步的了解。§1介绍现实对象和它的模型的关系，给出一些模型形式，说明什么是数学模型；§2阐述建立数学模型的重要意义；§3节通过几个示例说明用数学语言和数学方法表述和解决实际问题，即建立数学模型的过程；§4述建立数学模型的一般方法和步骤；§5介绍数学模型的特点及数学模型的分类；§6论建立数学模型能力的培养。
§1 从现实对象到数学模型人类生活在丰富多彩、变化万千的现实世界里，无时无刻不在运用智慧和力量去认识、利用、改造这个世界，从而不断地创造出日新月异、五彩缤纷的物质文明和精神文明。博览会常常是集中展示这些成果的场所之一，那些五光十色、精美绝伦的展品给我们留下了深刻的印象。工业博览会上，豪华、舒适的新型汽车叫人赞叹不已；农业博览会上，硕大、娇艳的各种水果令人流连忘返；科技展览厅里，大型水电站模型雄伟壮观，人造卫星模型高高耸立，清晰的数字和图表显示着电力工业的迅速发展，和整面墙壁一样大的地图上鲜明地标出了新建的铁路和新辟的航线，核电站工程的彩色巨照前，手持原子结构模型的讲解员深入浅出地介绍反应堆的运行机理；电影演播室里，播放着一部现代化炼钢厂实现生产自动控制的科技影片，其中既有火花四溅的钢坯浇铸情景，也有展示计算机管理和控制的框图、公式和程序。
参观博览会，像汽车、水果那些原封不动地从现实世界搬到展厅里的物品固然给人以亲切真实的感受，可是从开阔眼界、丰富知识的角度看，电站、卫星、铁路、钢厂……这些在现实世界被人们认识、构造、控制的对象，以它们的各种形式的模型——实物模型、照片、图表、公式、程序……汇集在人们面前，这些模型在短短几小时里所起的作用，恐怕是置身现实世界多少天也无法做到的。
与形形色色的模型相对应，它们在显示世界里的原始参照物通称为原型。本节先讨论原型和模型，特别是数学模型的关系，再介绍数学模型的意义。
原型和模型 原型(Prototupe)和模型(Model)是一对对偶体，原型指人们在现实世界里关心、研究或者从事生产、管理的实际对象。在科技领域通常使用系统(System)、过程(Process)等词汇、如机械系统、电力系统、生态系统、生命系统、社会经济系统，又如钢铁冶炼过程、导弹飞行过程、化学反应过程、污染扩散过程、生产销售过程、计划决策过程等。本书所述的现实对象、研究对象、实际问题等均指原型。模型则指为了某个特定目的将原型的某一部分信息简缩、提炼而构造的原型替代物。
这里特别强调构造模型的目的性。模型不是原型原封不动的复制品，原型有各个方面的各种层次的特征，而模型只要求反映与某种目的有关的那些方面和层次。一个原型，为了不同的目的可以有许多不同的模型。如放在展厅里的飞机模型应该在外形上逼真，但是不一定会飞。而参加航模竞赛的模型飞机要具有良好的飞行性能，在外观上不必苛求。至于在飞机设计、试制过程中用到的数学模型和计算机模拟，则只要在数量规律上真实反映飞机动态特性，毫不涉及飞机的实体。所以模型的基本特征是由构造模型的目的决定的。
我们已经看到模型有各种形式。用模型替代原型的方式来分类，模型可以分为物质模型（形象模型）和理想模型（抽象模型）。前者包括直观模型、物理模型等，后者包括思维模型、符号模型、数学模型等。
直观模型 指那些供展览用的实物模型，以及玩具、照片等，通常是把原型的尺寸按比例缩小或放大，主要追求外观上的逼真。这类模型的效果是一目了然的。
物理模型 主要指科技工作者为一定目的根据相似原理构造的模型，它不仅可以显示原型的外形或某些特征，而且可以用来进行模拟实验，间接地研究原型的某些规律。如波浪水箱中的舰艇模型用来模拟波浪冲击下舰艇的航行性能，风洞中的飞机模型用来试验飞机在气流中的空气动力学特性。有些现象直接用原型研究非常困难，更可借助于这类模型，如地震模拟装置、核爆炸反应模拟设备等。应注意验证原型与模型间的相似关系，以确定模拟实验结果的可靠性。物理模型常可得到实用上很有价值的结果，但也存在成本高、时间长、不灵活等缺点。
思维模型 指通过人们对原型的反复认识，将获取的知识以经验形式直接贮存于人脑中，从而可以根据思维或直觉作出相应的决策。如汽车司机对方向盘的操纵、一些技艺性较强的工种（如钳工）的操作，大体上是靠这类模型进行的。通常说的某些领导者凭借经验作决策也是如此。思维模型便于接受，也可以在一定条件下获得满意的结果，但是它往往带有模糊性、片面性、主观性、偶然性等缺点，难以对它的假设条件进行检验，并且不便于人们的互相沟通。
符号模型 是在一些约定或假设下借助于专门的符号、线条等，按一定形式组合起来描述原型。如地图、电路图、化学结构式等，具有简明、方便、目的性强及非量化等特点。
本书要专门讨论的数学模型则是由数字、字母或其它数学符号组成的，描述现实对象数量规律的数学公式、图形或算法。
什么是数学模型 其实你早在学习初等代数的时候就已经碰到过数学模型了。当然其中许多问题是老师为了教会学生知识而人为设置的。譬如你一定解过这样的所谓“航行问题”：
甲乙两地相距750 km，船从甲到乙顺水航行需30 h，从乙到甲逆水航行需50 h，问船速、水速各若干？
用，分别代表船速和水速，可以列出方程
,
实际上，这组方程就是上述航行问题的数学模型。列出方程，原问题已转化为纯粹的数学问题。方程的解km/h， km/h，最终给出了航行问题的答案。
当然，真正实际问题的数学模型通常要复杂得多，但是建立数学模型的基本内容已经包含在解这个代数应用题的过程中了。那就是：根据建立数学模型的目的和问题的背景做出必要的简化假设（航行中设船速和水速为常数）；用字母表示待求的未知量（，代表船速和水速）；利用相应的物理或其它规律（匀速运动的距离等于速度乘以时间），列出数学式子（二元一次方程组）；求出数学上的解答（，）；用这个答案解释原问题（船速和水速分别为20 km/h和5 km/h）；最后还要用实际现象来验证上述结果。
一般地说，数学模型可以描述为，对于现实世界的一个特定对象，为了一个特定目的，根据特有的内在规律，做出一些必要的简化假设，运用适当的数学工具，得到的一个数学结构。
需要指出，本书的重点不在于介绍现实对象的数学模型(Mathematical Model)是什么样子，而是要讨论建立数学模型(Mathematical Modelling)的全过程，建立数学模型下面简称为数学建模或建模。
与数学建模有密切关系的数学模拟，主要指运用数字式计算机的计算机模拟(Computer Simulation)。它根据实际系统或过程的特性，按照一定的数学规律用计算机程序语言模拟实际运行状况，并依据大量模拟结果对系统或过程进行定量分析。例如通过各种工件在不同机器上按一定工艺顺序加工的模拟，可以分析车辆在路段上的分布特别是堵塞的状况，与用物理模型的模拟实验相比，计算机模拟有明显的优点；成本低、时间短、重复性高、灵活性强。有人把计算机模拟作为建立数学模拟的手段之一，但是数学模拟在某种意义下描述了对象内在特性的数量关系，其结果容易推广，特别是得到了解析形式答案时，更容易推广。而计算机模拟则完全模仿对象的实际演变过程，难以从得到的数字结果分析对象的内在规律。当然，对于那些因内部机理过于复杂，目前尚难以建立数学模型的实际对象，用计算机模拟获得一些定量结果，可称是解决问题的有效手段。
§2 数学建模的重要意义数学，作为一门研究现实世界数量关系和空间形式的科学，在它产生和发展的历史长河中，一直是和人们生活的实际需要密切相关的。作为用数学方法解决实际问题的第一步，数学建模自然有着与数学同样悠久的历史。两千多年以前创立的欧几里德几何，17世纪发现的牛顿万有引力定律，都是科学发展史上数学建模的成功范例。
进入20世纪以来，随着数学以空间的广度和深度向一切领域的渗透，和电子计算机的出现与飞速发展，数学建模越来越受到人们的重视，可以从以下几方面来看数学建模在现实世界中的重要意义。
1）在一般工程技术领域，数学建模仍然大有用武之地在以声、光、热、力、电这些物理学科为基础的诸如机械、电机、土木、水利等工程技术领域中，数学建模的普遍性和重要性不言而喻。虽然这里的基本模型是已有的，但是由于新技术、新工艺的不段涌现，提出了许多需要用数学方法解决的新问题；高速、大型计算机的飞速发展，使得过去即便有了数学模型也无法求解的课题（如大型水坝的应力计算，中长期天气预报等）迎刃而解；建立在数学模型和计算机模拟基础上的CAD技术，以其快速、经济、方便等优势，大量地替代了传统工程设计中的现场实验、物理模拟等手段。
2）在高薪技术领域，数学建模几乎是必不可少的工具无论是发展通讯、航天、微电子、自动化等高新技术本身，还是将高新技术用于传统工业去创造新工艺、开发新产品，计算机技术支持下的建模和模拟都是经常使用的有效手段。数学建模、数值计算和计算机图形学等相结合形成的计算机软件，已经被固化于产品中，在许多高新技术领域起着核心作用，被认为是高新技术的特征之一。在这个意义上，数学不再仅仅作为一门科学，是许多技术的基础，而且直接走向了技术的前台。国际上一位学者就提出了“高技术本质上一种数学技术”的观点。
3）数学迅速进入一些新领域，为数学建模开拓了许多新的处女地随着数学向诸如经济、人口、生态、地质等所谓非物理领域的渗透，一些交叉学科如计量经济学、人口控制论、数学生态学、数学地质学等应运而生。这里一般地说不存在作为支配关系的物理定律，当用数学方法研究这些领域中的定量关系时，数学建模就成为首要的、关键的步骤和这些学科发展与应用的基础。在这些领域里建立不同类型、不同方法、不同深浅程度的模型的余地相当大，为数学建模提供了广阔的新天地。马克思说过：“一门科学只有成功地运用数学时，才算达到了完善的地步”。展望21世纪，数学必将大踏步地进入所有学科，数学建模将迎来蓬勃发展的新时期。
今天，在国民经济和社会活动的以下诸多方面，数学建模都有着非常具体的应用。
分析与设计 例如描述药物浓度在人体内的变化规律以分析药物的疗效；建立跨音速流和激波的数学模型，用数值模拟设计新的飞机翼型。
预报与决策 生产过程中产品质量指标的预报、气象预报、人口预报、经济增长预报等等，都要有预报模型；使经济效益最大的价格策略、使费用最少的设备维修方案，是决策模型的例子。
控制与优化 电力、化工生产过程的最优控制、零件设计中的参数优化，要以数学模型为前提。建立大系统控制与优化的数学建模，是迫切需要和十分棘手的课题。
规划与管理 生产计划、资源配置、运输网络规划、水库优化调度，以及排队策略、物资管理等，都可以用数学规划模型解决。
数学建模与计算机技术的关系密不可分。一方面，像新型飞机设计、石油勘探数据处理中的数学模型的求解当然离不开巨型计算机，而微型电脑的普及更使数学建模逐步进入人们的日常活动。比如当一位公司经理根据客户提出的产品数量、质量、交货期等要求，用手提电脑与客户进行价格谈判时，您不会怀疑他的电脑中贮存了由公司的各种资源、产品工艺流程及客户需求等数据研制的数学模型——快速报价系统和生产计划系统。另一方面，以数字化为特征的信息正以爆炸之势涌入计算机，去伪存真、归纳整理、分析现象、显示结果……，计算机需要人们给它以思维的能力，这些当然要求助于数学模型，所以把计算机技术与数学建模在知识经济中的作用比喻为如虎添翼，是恰如其分的。
美国科学院一位院士总结了将数学科学转化为生产力过程中的成功和失败，得出了“数学是一种关键的、普遍的、可以应用的技术”的结论，认为数学“由研究到工业领域的技术转化，对加强经济竞争力具有重要意义”，而“计算和建模重新成为中心课题，它们是数学科学技术转化的主要途径”。
§3.1 建模示例之一 椅子能在不平的地面上放稳吗本节和下面两节将给出三个数学建模的例子，重点说明如何作出合理的、简化的假设，用数学语言确切地描述实际问题，以及模型的结果怎样解释实际现象。
本节讨论的问题来源于日常生活中一件普通的事实；把椅子往不平的地面上一放，通常只有三只脚着地，放不稳，然而只需稍微挪动几次，就可以使四只脚同时着地，放稳了。这个看来似乎与数学无关的现象能用数学语言给予表述，并用数学工具来证实吗？让我们试试看[32].
模型假设 对椅子和地面都要作一些必要的假设：
椅子四条腿一样长，椅脚与地面接触可视为一个点，四脚的连线呈正方形。
地面高度是连续变化的，沿任何方向都不会出现间断（没有像台阶那样的情况），即地面可视为数学上的连续曲面。
对于椅脚的间距和椅脚的长度而言，地面是相对平坦的，使椅子在任何位置至少有三只脚同时着地。
假设1显然是合理的。假设2相当于给出了椅子能放稳的条件，因为如果地面高度不连续，譬如在有台阶的地方是无法使四只脚同时着地的。至于假设3是要排除这样的情况；地面上与椅脚间距和椅腿长度的尺寸大小相当的范围内，出现深沟或凸峰（即使是连续变化的），致使三只脚无法同时着地。
模型构成
中心问题是用数学语言把椅子四只脚同时着地的条件和结论表示出来。
首先用变量表示椅子的位置，由于椅脚的连线呈正方形，以中心为对称点，正方形绕中心的旋转正好代表了椅子的位置改变，于是可以用旋转角度这一变量来表示椅子的位置。
其次要把椅脚着地用数学符号表示出来，如果用某个变量表示椅脚与地面的竖直距离，当这个距离为0时，表示椅脚着地了。椅子要挪动位置说明这个距离是位置变量的函数。
由于正方形的中心对称性，只要设两个距离函数就行了，记A、C两脚与地面距离之和为，B、D两脚与地面距离之和为，显然、，由假设（2）知f、g都是连续函数，再由假设（3）知、至少有一个为0。当时，不妨设，这样改变椅子的位置使四只脚同时着地，就归结为如下命题：
命题 已知、是的连续函数，对任意，*=0，且，则存在，使。
模型求解将椅子旋转，对角线AC和BD互换，由可知。令，则，由f、g的连续性知h也是连续函数，由零点定理，必存在使，，由，所以。
评 注模型巧妙在于用已知变量表示椅子的位置，用的两个函数表示椅子四脚与地面的距离。利用正方形的中心对称性及旋转并不是本质的，同学们可以考虑四脚呈长方形的情形。
§3.2 建模示例之二 商人们怎样安全过河三名商人各带一个随从乘船渡河，一只小船只能容纳二人，由他们自己划行。随从们密约，在河的任一岸，一旦随从的人数比商人多，就杀人越货。但是如何乘船渡河的大权掌握在商人们的手中。商人们怎样才能安全渡河呢？
对于这类智力游戏经过一番逻辑思索是可以找出解决办法的。这里用数学模型求解，一是为了给出建模的示例，二是因为这类模型可以解决相当广泛的一类问题，比逻辑思索的结果容易推广。
由于这个虚拟的问题已经理想化了，所以不必再作假设。安全渡河问题可以视为一个多步决策过程。每一步，即船由此岸驶向彼岸或从彼岸驶回此岸，都要对船上的人员（商人、随从各几人）作出决策，在保证安全的前提下（两岸的随从数都不比商人多），在有限步内使全部人员过河。用状态（变量）表示某一岸的人员状况，决策（变量）表示船上的人员状况，可以找出状态随决策变化的规律。问题转化为在状态的允许变化范围内（即安全渡河条件），确定每一步的决策，达到渡河的目标。
模型构成 记第次渡河前此岸的商人数为，随从数为，，,将二维向量定义为状态。安全渡河条件下的状态集合称为允许状态集合，记作.
 (1)
不难验证，对此岸和彼岸都是安全的。
记第次渡船上的商人数为，随从数为,将二维向量定义为决策。允许决策集合记作，由小船的容量可知
 (2)
因为为奇数时船从此岸驶向彼岸，为偶数时船由彼岸驶回此岸，所以状态随决策变化的规律是
 (3)
(3)式称状态转移律。这样，制定安全渡河方案归结为如下的多步决策模型：
求决策，使状态按照转移率(3)，由初始状态经有限步到达状态.
模型求解 根据(1)~ (3)式编一段程序用计算机求解上述多步决策问题是可行的。不过对于商人和随从人数不大的简单状况，用图解法这个模型更为方便。

图2 安全渡河问题的图解法在平面坐标系上画出图2那样的方格，方格点表示状态,允许状态集合是用圆点标出的10个格子点。允许决策是沿方格线移动1或2格，为奇数时向左、下方移动，为偶数是向右、上方移动。要确定一系列的使由经过那些圆点最终移至原点.
图2给出了一种移动方案，经过决策，最终有,这个结果很容易翻译成渡河的方案。
评注 这里介绍的是一个规格化的方法，所建立的多步决策模型可以用计算机求解，从而具有推广的意义。譬如当商人和随从人数增加或小船的容量加大时，靠逻辑思考就困难了，而用这种模型则仍可方便地求解。读者不妨考虑四名商人各带一个随从的情况（小船同前）。
适当地设置状态和决策，确定状态转移律，建立多步决策模型，是有效地解决很广泛的一类问题的方法，在以后的章节中还要用到。
§3.3 建模示例之三 如何预报人口的增长人类社会进入20世纪以来，在科学技术和生产力飞速发展的同时，世界人口也以空前的规模增长。统计数据显示：
年
1625
1830
1930
1960
1974
1987
1999
人口（亿）
5
10
20
30
40
50
60
可以看出，人口每增加十亿的时间，由一百年缩短为十二三年。我们赖以生存的地球，已经携带着它的60亿子民踏入21世纪。
长期以来，人类的繁殖一直在自发地进行着。只是由于人口数量的迅速膨胀和环境质量的急剧恶化，人们才猛然醒悟，开始研究人类和自然的关系、人口数量的变化规律，以及如何进行人口控制等问题。
我国是世界第一人口大国，地球上每九个人中就有一个中国人。在20世纪的一段时间内我国人口的增长速度过快，请看：
年
1908
1933
1953
1964
1982
1990
2000
人口（亿）
3.0
4.7
6.0
7.2
10.3
11.3
12.95
有效地控制我国人口的增长，不仅是使我国全面进入小康社会、到21世纪中叶建成富强民主文明的社会主义国家的需要，而且对于全人类社会的美好理想来说，也是我们义不容辞的责任。
认识人口数量的变化规律，建立人口模型，做出准确的预报，是有效控制人口增长的前提。长期以来人们在这方面作了不少工作，下面介绍两个最基本的人口模型，并利用表1给出的近两个世纪的美国人口统计数据（以百万为单位），对模型作检验，最后用它预报2010年美国的人口。
年人口
1790
3.9
1800
5.3
1810
7.2
1820
9.6
1830
12.9
1840
17.1
1850
23.2
1860
31.4
年人口
1870
38.6
1880
50.2
1890
62.9
1900
76.0
1910
92.0
1920
106.5
1930
123.2
1940
131.7
年人口
1950
150.7
1960
179.3
1970
204.0
1980
226.5
1990
251.4
2000
281.4
表1 美国人口统计数据
1) 指数增长模型最简单的人口增长模型使人所共识的：记今年人口为，年后人口为，年增长率为，则
 (1)
显然，这个公式的基本条件是年增长率保持不变。
二百多年前英国人口学家马尔塞斯（Malths，1766—1834）调查了英国一百多年的人口统计资料，得出了人口增长率不变的假设，并据此建立了著名的人口指数增长模型。
模型建立 记时刻的人口为，当考察一个国家或一个较大地区的人口时，是一个很大的整数。为了利用微积分这一数学工具，将视为连续、可微函数。记初始时刻（）的人口为,假设人口增长率为常数，即单位时间内的增量等于乘以,考虑到时间内人口的增量，显然有

令，得到满足微分方程：
， (2)
有这个方程很容易解出
 (3)
时(3)式表示人口将按指数规律随时间无限增长，称为指数增长模型。
参数估计 (3)式的参数和可以用表1的数据估计。为了利用简单的线性最小二乘法，将(3)式取对数，可得
，， (4)
以1790年至1900年的数据拟合(4)式，用MATLAB软件计算可得=0.2743/10年，=4.1884,以全部数据（1790年至2000年）拟合(4)式，得=0.2022/10年，=6.0450.
结果分析 用上面得到的参数和代入(3)式，将计算结果与实际数据作比较。表2中计算人口是用1790年至1900年的数据拟合的结果，计算人口是用全部数据（1790年至2000年）拟合的结果，图3、图4是它们的图形表示（+号是实际数据，曲线是计算结果）。可以看出，用这个模型基本上能够描述十九世纪以前美国人口的增长，但是进入20世纪后，美国人口增长明显变慢，这个模型就不合适了。
年
实际人口
计算人口
计算人口
1790
3.9
4.2
6.0
1800
5.3
5.5
7.4
1810
7.2
7.2
9.1
1820
9.6
9.5
11.1
1830
12.9
12.5
13.6
1840
17.1
16.5
16.60
1850
23.2
21.7
20.30
1860
31.4
28.6
24.90
1870
38.6
37.6
30.5
1880
50.2
49.5
37.3
1890
62.9
65.1
45.7
1900
76.0
85.6
55.9
1910
92.0
68.4
1920
106.5
83.7
1930
123.2
102.5
1940
131.7
125.5
1950
150.7
153.6
1960
179.3
188.0
1970
204.0
230.1
1980
226.5
281.7
1990
251.4
344.8
2000
281.4
422.1
表2 指数增长模型拟合美国人口数据的结果

图3 指数增长模型拟合图形

图4 指数增长模型拟合图形练习 用1900至2000年的数据拟合指数增长模型，计算并作图，观察结果。
历史上，指数增长模型与十九世纪以前欧洲一些地区人口统计数据可以很好地吻合，迁往加拿大的欧洲移民后代人口也大致符合这个模型。另外，用它做短期人口预测可以得到较好的结果。显然，这是因为在这些情况下，模型的基本设计——人口增长率是常数——大致成立。
但是长期以来，任何地区的人口都不可能无限增长，即指数模型不能描述、也不能预测较长时期的人口演变过程。这是因为，人口增长率事实上是在不断地变化着。排除灾难、战争等特殊时期，一般说来，当人口较少时，增长较快，即增长率较大；人口增加到一定数量以后，增长就会慢下来，即增长率变小。表3是用数值微分的三点公式计算的美国人口增长率（﹪/年），可以看到，进入20世纪后增长率明显下降。用平均增长率作为，用指数增长模型描述美国人口的变化，结果当然会与表1的统计数据相差很大。
年增长率
1790
2.95
1800
3.11
1810
2.99
1820
2.97
1830
2.91
1840
3.01
1850
3.08
1860
2.45
年增长率
1870
2.44
1880
2.42
1890
2.05
1900
1.91
1910
1.66
1920
1.46
1930
1.02
1940
1.04
年增长率
1950
1.58
1960
1.49
1970
1.16
1980
1.05
1990
1.09
2000
1.16
表3 美国人口增长率（﹪/年）
看来，为了使人口预报特别是长期预报更好地符合实际情况，必须修改指数增长模型关于人口增长率是常数这个基本假设[21].
2) 阻滞增长模型（Logistic模型）
模型建立 分析人口增长到一定数量后增长率下降的主要原因，人们注意到，自然资源、环境条件等因素对人口的增长起着阻滞作用，并且随着人口的增长，阻滞作用越来越大。所谓阻滞增长模型就是考虑到这个因素，对指数增长模型的基本假设进行修改后得到的。
阻滞作用体现在对人口增长率的影响上，使得随着人口数量的增加而下降。若将表示为的函数，则它应是减函数。于是方程(2)写作
， (5)
对的一个最简单的假定是，设为的线性函数，即
 (6)
这里称固有增长率，表示人口很少时（理论上是=0）的增长率。为了确定系数的意义，引入自然资源和环境条件所能容纳的最大人口数量，称人口容量。当时人口不再增长，即增长率应，代入(6)式得，于是(6)式为
 (7)
(7)式的另一种解释是，增长率与人口尚未实现部分的比例成正比，比例系数为固有增长率.
将(7)代入方程(5)得
， (8)
方程(8)右端的因子体现人口自身的增长趋势，因子则体现了资源和环境对人口增长的阻滞作用。显然，越大，前一因子越大，后一因子越小，人口增长是两个因子共同作用的结果。
如果以为横轴，为纵轴作出方程(8)的图形（图5），可以分析人口增长速度随着的增加而变化的情况，从而大致地看出的变化规律。
练习 根据图5 与的关系，分析随的变化情况：较小（从而较小）时和较大（从而较大）时的增长速度有何不同，多大时人口增长最快，时？等，由此你能大致画出的图形吗。
实际上，方程(8)可以用分离变量法求解得到
 (9)
读者可以用计算机画出(9)式的图形，它是一条S形曲线（图6），增加得先快后慢，时，拐点在,你在上面的联系中画的图形与这个图形一样吗。
参数估计 为了利用简单的线性最小二乘法估计这个模型的参数和，我们不用(9)式，而将方程(8)表为
 (10)
上式左端可以从表1的数据用数值微分计算，右端对参数，是线性的。我们利用1860年至1990年的数据（去掉个别异常数据），用MATLAB软件计算得到=0.2557/10年，=392.0886.
参数估计也可以借助专家的经验。例如，某些人口学家估计世界人口的固有增长率=0.029，又知道世界人口在1960年为29.8亿时，增长率是1.85﹪，即=0.0185，于是按照方程(8)，世界人口容量为=29.8(1-0.0185/0.029)=82.3 亿。实际上，20世纪70年代世界人口为40亿左右时增长率达到最大，然后开始下降。注意到阻滞增长模型中时最大，可以看出上述结果的一致性。
结果分析 用上面得到的参数和代入(9)式，将计算结果与实际数据作比较，得表4和图7
年
实际人口
计算人口
1790
3.9
3.9
1800
5.3
5.0
1810
7.2
6.5
1820
9.6
8.3
1830
12.9
10.7
1840
17.1
13.7
1850
23.2
17.5
1860
31.4
22.3
1870
38.6
28.3
1880
50.2
35.8
1890
62.9
45.0
1900
76.0
56.2
1910
92.0
69.7
1920
106.5
85.5
1930
123.2
103.9
1940
131.7
124.5
1950
150.7
147.2
1960
179.3
171.3
1970
204.0
196.2
1980
226.5
221.2
1990
251.4
245.3
表4 阻滞增长模型拟合美国人口数据的结果

图7 阻滞增长模型拟合图形（以1790年为起点）
可以看出，用这个模型拟合时虽然中间一段（19世纪中叶到20世纪中叶）不大好，但是最后一段（二十世纪中叶以后）吻合的不错。
模型检验 在估计阻滞增长模型的参数时没有用2000年的实际数据，是为了用它作模型检验。我们用模型计算2000年的人口，与已知的实际数据（281.4百万）比较，来检验模型是否合适。
为简单起见，可利用和方程(7)作如下计算

得到=274.5百万，与实际数据的误差约2.5﹪，可以分为该模型是相当满意的。
人口预报 应将2000年的实际数据加进去重新估计参数，可得=0.2490/10年，=433.9886,然后再利用模型检验中的计算方法预报美国2010年的人口，得到年的人口，得到=306.0百万。这个预报结果的准确性如何。让我们拭目以待。
由方程(8)表示的阻滞增长模型，是荷兰生物学家Verhulst19世纪中叶提出的。它不仅能够大体上描述人口及许多物种数量（如森林中的树木、鱼塘中的鱼群等）的变化，而且在社会经济领域也有广泛的应用，例如耐用消费平的销售就可以用它来描述。基于这个模型能够描述一些事物符合逻辑的客观规律，人们常称它为Logistic模型，本书以后章节将多次用到它[31].
评注 用数学工具描述人口变化规律，关键是对人口增长率作出合理、简化的假定。阻滞增长模型就是将指数增长模型关于人口增长率是常数的假设进行修正后得到的。可以想到，影响增长率的出生率与死亡率与年龄有关，所以，更合乎实际的人口模型应该考虑年龄的因素。后面的讨论这样的模型。
参数估计和模型检验是建模的重要步骤，本节用到的线性最小二乘法是参数估计（如果是基于数据的，又称数据拟合）的基本方法，读者应该知道它的原理，并掌握它的软件实现。
§4 数学建模的基本方法和步骤数学建模面临的实际问题是多种多样的，建模的目的不同、分析的方法不同、采用的数学工具不同，所得模型的类型也不同，我们不能指望归纳出若干条准则，适用于一切实际问题的数学建模方法。下面所谓基本方法不是针对具体问题而是从方法论的意义上讲的。
数学建模的基本方法一般说来建模方法大体上可分为机理分析和测试分析两种。机理分析是根据对客观事物特性的认识，找出反映内部机理的数量规律，建立的模型常有明确的物理或现实意义。前面几个示例都是用的机理分析。测试分析将研究对象看做一个“黑箱”系统（意义是它的内部机理看不清楚），通过对系统输入、输出数据的测量和统计分析，按照一定的准则找出与数据拟合得最好的模型。
面对一个实际问题用哪一种方法建模，主要取决与人们对研究对象的了解程度和建模目的。如果掌握了一些内部机理的知识，模型也要求具有反映内在特征的物理意义，建模就应以机理分析为主。而如果对象的内部规律基本上不清楚，模型也不需要反映内部特性（例如仅用于对输出作预报），那么就可以用测试分析。
对于许多实际问题还常常将两种方法结合起来建模，即用机理分析建立模型的结构，用测试分析确定模型的参数。1.5节建立的人口模型就是这种情况。
机理分析当然要针对具体问题来做，不可能有统一的方法。因而主要是通过实例研究(Case studies)来学习。测试分析有一套完整的数学方法，第10章统计回归模型是其中的一小部分。以动态系统为主的测试分析称为系统辨识(System Identification)，是一门专门学科。本书以后所说的数学建模主要指机理分析。
数学建模的一般步骤建模要经过哪些步骤并没有一定的模式，通常与问题性质、建模目的等有关。下面介绍的是机理分析方法建模的一般过程，如图8所示

图8 数学建模步骤示意图模型准备 了解问题的实际背景，明确建模目的，搜集必要的信息如现象、数据等，尽量弄清楚对象的主要特征，形成一个比较清晰的“问题”，由此初步确定用哪一类模型。情况明才能方法对。在模型准备阶段要深入调查研究，虚心向实际工作者请教，尽量掌握第一手资料。
模型假设 根据对象的特征和建模目的，抓住问题的本质，忽略次要因素，作出必要的、合理的简化假设。对于建模的成败这是非常重要和困难的一步。假设作得不合理或太简单，会导致错误的或无用的模型；假设作得过分详细，试图把复杂对象的众多因素都考虑进去，会使你很难或无法继续下一步的工作，常常需要在合理与简化之间作出恰当的折衷。通常，作假设的依据，一是处于对问题内在规律的认识，二是来自对现象、数据的分析，以及二者的综合。想象力、洞察力、判断力，以及经验，在模型假设中起着重要的作用。
模型构成 根据所作的假设，用数学的语言、符号描述对象的内在规律，建立包含常量、变量等的数学模型，如优化模型、微分方程模型、差分方程模型、跳的模型等。这里除了需要一些相关学科的专门知识外，还常常需要较为广阔的应用数学方面的知识。要善于发挥想象力，注意使用类比法，分析对象与熟悉的其它对象的共性，借用已有的模型。建模时还应遵循的一个原则是：尽量采用简单的数学工具，因为你的模型总是希望更多的人了解和使用，而不是只供少数专家欣赏。
模型求解 可以采用解方程、画图形、优化方法、数值计算、统计分析等各种数学方法，特别是数学软件和计算机软件。
模型分析 对求解结果进行数学上的分析，如结果的误差分析、统计分析、模型对数据的灵敏性分析、对假设的强健性分析等。
模型检验 把求解和分析结果翻译回到实际问题，与实际的现象、数据比较，检验模型的合理性与适用性。如果结果与实际不符，问题常常出在模型假设上，应该修改、补充假设，重新建模，如图8中的虚线所示。这一步对于模型是否真的有用非常关键，要以严肃认真的态度对待。有些模型要经过几次反复，不断完善，直到检验结果获得某种程度上的满意。
模型应用 应用的方式与问题性质建模目的及最终的结果有关，一般不属于本书讨论的范围。
应当指出，并不是所有问题的建模都要经过这些步骤，有时各步骤之间的界限也不那么分明，建模时不要拘泥于形式上的按部就班，本书的实例就采取了灵活的表达形式。
数学建模的全过程
从前面的几个建模示例以及一般步骤的分析，可以将数学建模的过程分为表述、求解、解释、验证几个阶段，并且通过这些阶段完成从现实对象到数学模型，再从数学模型回到现实对象的循环，如图9所示。

图9 数学建模的全过程表述是将现实问题“翻译”成抽象的数学问题，属于归纳法。数学模型的求解则属于演绎法。归纳是依据个别现象推出一般规律；演绎是按照普遍原理考察特定对象，导出结论。因为任何事物的本质都要通过现象来反映，必然要透过偶然来表露，所以正确的归纳不是主观、盲目的，而是有客观基础的，但也往往是不精细的、带感性的，不宜直接检验其正确性。演绎利用严格的逻辑推理，对解释想象、作出科学预见具有重要意义，但是它要以归纳的结论作为公理化形式的前提，只能在这个前提下保证其正确性。因此，归纳和演绎是辩证统一的过程：归纳是演绎的基础，演绎是归纳的指导[51]。
解释是把数学模型的解答“翻译”回到现实对象，给出分析、预报、决策或者控制的结果、最后，作为这个过程的重要一环，这些结果需要用实际的信息加以验证。
图9也揭示了现实对象和数学模型的关系。一方面，数学模型是将现象加以归纳、抽象的产物，它源于现实，又高于现实。另一方面，只有当数学建模的结果经受住现实对象的检验时，才可以用来指导实际，完成实践——理论——实践这一循环。
§5 数学建模的特点和分类数学建模是利用数学工具解决实际问题的重要手段，得到的模型有许多有点，也有一些弱点。下面归纳出数学模型的若干特点，以期读者在学习过程中逐步领会[26]。
数学模型的特点模型的逼真性和可行性 一般说来总是希望模型尽可能逼近研究对象，但是一个非常逼真的模型在数学上常常是难于处理的，因而不容易达到通过建模对现实对象进行分析、预报、决策或者控制的目的，即实用上不可行。另一方面，越逼真的模型常常越复杂，即使数学上能处理，这样的模型应用时所需要的“费用”也相当高，而“费用”不一定与复杂模型取得的“效益”相匹配。所以建模时往往需要在建模的逼真性与可行性，“费用”与“效益”之间做出折中和抉择。
模型的渐进性 稍微复杂一些的实际问题的建模通常不可能一次成功，要经过上一节描述的建模过程的反复迭代，包括由简到繁，也包括删繁就简，以获得越来越满意的模型。在科学发展过程中随着人们认识和实践能力的提高，各门学科中的数学模型也存在着一个不断完善或者推陈出新的过程。从19世纪力学、热学、电学等许多学科由牛顿力学的模型主宰，到20世纪爱因斯坦相对论模型的建立，是模型渐进性的明显例证。
模型的强健性 模型的结构和参数常常是由模型假设及对象的信息如观测数据确定的，而假设不可能太准确，观测数据也是允许有误差的。一个好的模型应该具有下述意义的强健性：当模型假设改变时，可以导出模型结构的相应变化；当观测数据有微小改变时，模型参数也只有相应的微小变化。
模型的可转移性 模型是现实对象抽象化、理想化的产物，它不为对象的所属领域所独有，可以转移到另外的领域。在生态、经济、社会等领域内建模就常常借用物理领域中的模型。模型的这种性质显示了它的应用的极端广泛性。
模型的非预制性 虽然已经发展了许多应用广泛的模型，但是实际问题是各种各样、变化万千的，不可能要求把各种模型做成预制品供你再建模时使用。模型的折中非预制性使得建模本身常常是事先没有答案的问题（Open-end problem）。在建立新的模型的过程中甚至会伴随着新的数学方法或数学概念的产生。
模型的条理性 从建模的角度考虑问题可以促使人们对现实对象的分析更全面、更深入、更具条理性，这样即使建立的模型由于种种原因尚未达到实用程度，对问题的研究也是有力的。
模型的技艺性 建模的方法与其他一些数学方法如方程解法、规划问题解法等是根本不同的，无法归纳出若干条普遍适用的建模准则和技巧。有人说，建模目前与其说是一门技术，不如说是一种艺术，是技艺性很强的技巧。经验、想象力、洞察力、判断力以及直觉、灵感等在建模过程中起得作用往往比一些具体的数学知识更大。
模型的局限性 这里有几方面的含义。第一，由数学模型得到的结论虽然具有通用性和精确性，但是因为模型是现实对象简化、理想化的产物，所以一旦将模型的结论应用于实际问题，就回到了现实世界，那些被忽视、简化的因素必须考虑，于是结论的通用性和精确性只是相对的和近似的。第二，由于人们认识能力和科学技术包括数学本身发展水平的限制，还有不少实际问题很难得到有着实用价值的数学模型。如一些内部机理复杂、影响因素众多、测量手段不够完善、技艺性较强的生产过程，像生铁冶炼过程，常常需要开发专家系统，与建立数学模型相结合才能获得较满意的应用结果。专家系统是一种计算机软件系统，它总结专家的知识和经验，模拟人类的逻辑思维过程，建立了若干规则和推理途径，主要是定性地分析各种实际现象并作出判断。专家系统可以看成计算机模拟的新发展。第三，还有些领域中的问题今天尚未发展到用建模方法寻求数量规律的阶段，如中医诊断过程，目前所谓计算机辅助诊断也是属于总结著名中医的丰富临床经验的专家系统。
数学模型的分类数学模型可以按照不同的方式分类，下面介绍常用的几种。
1,按照模型的应用领域（或所属学科）分。如人口模型、交通模型、环境模型、生态模型、城镇模型、城镇规划模型、水资源模型、再生资源利用模型、污染模型等。范畴更大一些则形成许多边缘学科如生物数学、医学数学、地质数学、数量经济学、数学社会学等。
2,按照建立模型的数学方法（或所属数学分支）分。如初等模型、几何模型、微分方程模型、统计回归模型、数学规划模型等。
3,按照模型的表现特性又有几种分法：
确定性模型和随机性模型 取决于是否考虑随机因素的影响。近年来随着数学的发展，又有所谓突变性模型和模糊性模型。
静态模型和动态模型 取决于是否考虑时间因素引起的变化。
线性模型和非线性模型 取决于模型的基本关系，如微分方程是否是线性的离散模型和连续模型 模型中的变量（主要是时间变量）取为离散还是连续的。
虽然从本质上讲大多数实际问题是随机性的、动态的、非线性的，但是由于确定性、静态、线性模型容易处理，并且往往可以作为初步的近似来解决问题，所以建模时常先考虑确定性、静态、线性模型。连续模型便于利用微积分方法求解析解，作为理论分析，而离散模型便于在计算机上作数值计算，所以用哪种模型要看具体问题而定。在具体的建模过程中将连续模型离散化，或将连续变量视作连续的，也是常采用的方法。
4,按照建模目的分。有描述模型、预报模型、优化模型、决策模型、控制模型等。
5,按照对模型结构的了解程度分。有所谓白箱模型、灰箱模型、黑箱模型。这是把研究对象比喻成一只箱子里的机关，要通过建模来解释它的奥妙。白箱主要包括用力学、热学、电学等一些机理相当清楚的学科描述的现象以及相应的工程技术问题，这方面的模型大多已经基本确定，还需深入研究的主要是优化设计和控制等问题了。灰箱主要指生态、气象、经济、交通等领域中机理尚不十分清楚的现象，在建立和改善模型方面都还不同程度地有许多工作要做。至于黑箱则主要指生命科学和社会科学等领域中一些机理（数量关系方面）很不清楚的现象。有些工程技术问题虽然主要基于物理、化学原理，但由于因素众多、关系复杂和观测困难等原因也常作为灰箱或黑箱模型处理。当然，白、灰、黑之间并没有明显的界限，而且随着科学技术的发展，箱子的“颜色”必然是逐渐由暗变亮的。
§6 数学建模能力的培养在详细分析了建立数学模型的全过程和数学模型的特点以后，我们看到用建模方法解决实际问题，首先是用数学语言表述问题及构造模型，其次才是用数学工具求解构成的模型。绝大多数数学课程如微积分、线性代数、概率论、计算方法等都是讲授某一专门知识和培养数学运算、逻辑推理能力的，这些数学技巧主要用来求解数学模型。用数学语言表述问题，包括模型假设、模型构造等，除了要有广博的知识（包括数学知识和各种实际知识）和足够的经验之外，特别需要丰富的想象力和明锐的洞察力。
想象力指人们在原有知识的基础上，将新感知的形象与记忆中的形象相互比较、重新组合、加工处理，创造出新的形象，是一种形象思维活动。洞察力指人们在充分占有资源的基础上，经过初步分析能迅速抓住主要矛盾，舍弃次要因素，简化问题的层次，对可以用哪些方法解决面临的问题，以及不同方法的优劣做出判断[51]。
类比方法和理想化方法是建模中常用的方法，它们运用与想象力、洞察力有密切关系。类比法注意到研究对象与已熟悉的另一对象具有某些共性，比较二者相似之处已获得对研究对象的新认识。喧杂什么对象进行类比，比较哪些相似的属性，在一定程度上是靠想象进行的。将交通流与水流类比来建立交通流模型是这方面例子。理想化方法是从观察和经验中通过想象和逻辑思维，把对象简化、纯化，使其升华到理想状态，以期更本质地揭示对象的固有规律。在一定条件下把物体看作质点，把实际位置看做数学上的点、线等都是理想化的结果。
建模过程是一种创造性思维过程，除了想象、洞察、判断这些属于思维、逻辑思维范畴的能力之外，直觉和灵感往往也起着不可忽视的作用。直觉是人们对新事物本质的极敏锐的领悟、理解或推断，灵感指在人们有意识或下意识思考过程中迸发出来的猜测、思路或判断。二者都具有突发性，且思维者本人往往说不清它的来路和道理。当由于各种限制利用已有知识难以对研究对象作出有效的推理和判断时，凭借相似、类比、猜测，外推等思维方式及不完整、不连续、不严密的，带启发性的直觉和灵感，去“战略性”地认识对象，是人类创造性思维的特点之一，也是人脑比按程序逻辑工作的计算机、机器人的高明之处。历史上不乏在科学家的直觉和灵感的火花中诞生的假说、论证和定律。当然，直觉和灵感不是凭空产生的，它要求人们具有丰富的背景知识，对问题进行反复的思考和艰苦探索，对各种思维方法运用娴熟。相互讨论和思想交锋，特别是不同专业的成员之间的探讨，是激发直觉和灵感的重要因素。所以由各种专门人才组成的所谓团队工作方式（Team work）越来越受到重视。
前面说过，建模可以看成一门艺术。艺术在某种意义下是无法归纳出几条准则或方法的。一名出色的艺术家需要大量的观摩和前辈的指教，更需要亲身的实践。类似的，掌握建模这门艺术，培养想象力和洞察力，不外乎认真作好这样两条：第一，学习、分析、评价、改造别人作过的模型。首先要弄懂它，分析为什么这么做，然后找出它的优缺点，并尝试改进的方法。第二，要亲自动手，踏实地做几个实际题目。为了这个目的本书采用实例研究方法，一方面给出在各个应用领域不同数学方法建模的大量实例，另一方面通过习题提供若干题目让读者自己练习。实例研究方法虽然不能按照严密的逻辑结构去讨论问题，不能划定这些方法的实用范围，其得到的结果也并非无可置疑，但它却是我们学习建模以解决实际问题的一种很生动，有效的方法。
附：近几年全国大学生数学建模竞赛题
1997
A
零件的参数设计
B
最优截断切割问题
1998
A
投资的收益与风险
B
灾情巡视路线
1999
A
自动化车床管理
B
钻井布局
2000
A
DNA序列分析
B
钢管订购与运输
2001
A
血管的三维重建
B
公交车调度
2002
A
车灯线光源的优化设计
B
彩票中的数学
2003
A
SARS的传播
B
露天矿生产的车辆安排
2004
A
奥运会临时超市网点设计
B
电力市场的输电阻塞管理
2005
A
长江水质的评价和预测
B
关于DVD的在线租赁
2006
A
B
2007
A
B
课堂练习：把四只脚连线呈长方形的椅子往不平的地面上一放，通常只有三只脚着地，放不稳，然而有人认为只要稍挪动几次，就可以四脚着地，放稳了，对吗？（分析、建立模型并给出结论）
班级,姓名,学号,成绩,
第二篇 应用数学软件-MATLAB入门
在科学和工程应用中，往往要进行大量的科学计算，其中包括以矩阵为基础的数学计算；这些计算一般来说难以用手工精确和快捷地进行，而且众多工程问题一般只要求得到满足精度的近似解即可，从而借助于计算机编写相应的程序进行近似计算就显得很有必要。目前用Basic、Fortran和C编制计算程序较多，但其既需要对有关算法有深刻的了解，还需要熟练掌握所用语言的语法和编程技巧；这对较多科学和工程技术人员而言，同时具备这两方面的技能就很有难度；而且用上述语言编制程序不但复杂，一般需要大量的人力和物力，而且影响工作进程和效率，为此，美国Mathwork公司于1967年推出了“Matrix Laboratory”（即矩阵实验室，缩写为Matlab）软件包，并不断进行更新和扩充，目前已成为全球应用最广泛最流行的软件之一。
目前最新的7.1版本（windows环境）是一种功能强、效率高便于进行科学和工程计算的交互式软件包。其中包括：一般数值分析、矩阵运算、概率统计、建模与系统控制和优化等应用程序。集应用程序和图形于一体便于使用的集成环境中，在此环境下所解问题的Matlab语言表述形式和其数学表达形式相同，不需要按传统的方法编程，就可解决工程、科学计算和数学学科中的许多问题。不过，Matlab作为一种新的计算机语言，要想运用自如，充分发挥它的威力，必须先系统地学习它。首先应该相信的是，由于使用Matlab编程运算与人进行科学计算的思路和表达方式完全一致，从而学习Matlab语言不象学习其它高级语言——如Basic、Fortran和C等那样难以掌握。
§1 MATLAB语言基础
MATLAB建立在向量、数组和矩阵的基础上，使用方便，人机界面直观，输出结果可视化，矩阵是MATLAB的核心。
§1.1 变量与函数
§1.1.1 变量
MATLAB中变量的命名规则是：
（1）变量名必须是不含空格的单个词；
（2）变量名区分大小写；
（3）变量名最多不超过19个字符；
（4）变量名必须以字母打头，之后可以是任意字母、数字或下划线，变量名中不允许使用标点符号。
§1.1.2 特殊变量表
特殊变量
取 值
ans
用于结果的缺省变量名
pi
圆周率
eps
计算机的最小数，当和1相加就产生一个比1大的数
inf
无穷大，如1/0
NaN
不定量，如0/0
i，j
i=j=
realmin
最小可用正实数
realmax
最大可用正实数
§1.1.3 数学运算符号及标点符号
+
加法运算，适用于两个数或两个同阶矩阵相加
—
减法运算
*
乘法运算
.*
点乘运算
/
除法运算
./
点除运算
^
乘幂运算
.^
点乘幂运算
\
反斜杠表示左除
（1）MATLAB的每条命令后，若为逗号或无标点符号，则显示该条命令运行的结果；若命令后为分号，则运行结果不显示；
（2）“%” 后面所有文字为注释；
（3）“...”表示续行。
§1.1.4 数学函数
函 数
名 称
函 数
名 称
sin(x)
正弦函数
asin(x)
反正弦函数
cos(x)
余弦函数
acos(x)
反余弦函数
tan(x)
正切函数
atan(x)
反正切函数
abs(x)
绝对值
max(x)
最大值
min(x)
最小值
sum(x)
元素的总和
sqrt(x)
开平方
exp(x)
以e为底的指数
log(x)
自然对数
log10(x)
以10为底的对数
sign(x)
符号函数
fix(x)
取整
§1.1.5 M文件
MATLAB的内部函数是有限的，有时为了研究某一个函数的各种性态，需要为MATLAB定义新函数，为此必须编写函数文件,函数文件是文件名后缀为M的文件，这类文件的第一行必须是一特殊字符function开始;
格式为：function 因变量名=函数名（自变量名）
函数值的获得必须通过具体的运算实现，并赋给因变量,
M文件建立方法：（1） 在Matlab中，点:File->New->M-file
（2） 在编辑窗口中输入程序内容
（3） 点：File->Save，存盘，M文件名必须与函数名一致。Matlab的应用程序也以M文件保存。
例1.1 定义函数
（1）建立M文件：fun.m
function f=fun(x)
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2
（2）可以直接使用函数fun.m
例1.2 计算f(1,2),只需在Matlab命令窗口键入命令：
x=[1 2]
fun(x)
§1.2 数组
§1.2.1 创建简单的数组
x=[a b c d e f ] 创建包含指定元素的行向量
x=first：last 创建从first开始，加1计数，到last结束的行向量
x=first：increment：last 创建从first开始，加increment计数，last结束的行向量
x=linspace(first，last，n) 创建从first开始，到last结束，有n个元素的行向量
x=logspace(first，last，n) 创建从first开始，到last结束，有n个元素的以10为底的对数分隔行向量.
程序：
x=[1 2 3 4 5 8 7 18]
y=1:7
z=3:2:9
v=[y z]
u=linspace(2,9,11)
结果：
x =1 2 3 4 5 8 7 18
y =1 2 3 4 5 6 7
z =3 5 7 9
v =1 2 3 4 5 6 7 3 5 7 9
u =2.00 2.70 3.40 4.10 4.80 5.50 6.20 6.90 7.60 8.30 9.00
§1.2.2 数组元素的访问
（1）访问一个元素,x(i)表示访问数组x的第i个元素,
（2）访问一块元素,x(a,b,c)表示访问数组x的从第a个元素开始，以步长为b到第c个元素（但不超过c），b可以为负数，b缺省时为1,
（3）直接使用元素编址序号,x([a b c d]) 表示提取数组x的第a、b、c、d个元素构成一个新的数组[x(a) x(b) x(c) x(d)],
程序,
x=1:9
y=x(2:2:8)
z=[x(1) x(6) x(8)]
结果：
x=1 2 3 4 5 6 7 8 9
y=2 4 6 8
z=1 6 8
§1.2.3 数组的方向
前面例子中的数组都是一行数列，是行方向分布的,称之为行向量,数组也可以是列向量，它的数组操作和运算与行向量是一样的，唯一的区别是结果以列形式显示；
产生列向量有两种方法：
直接产生 例 c=[1；2；3；4]
转置产生 例 b=[1 2 3 4]; c=b’
说明：以空格或逗号分隔的元素指定的是不同列的元素，而以分号分隔的元素指定了不同行的元素。
§1.2.4 数组的运算
标量-数组运算数组对标量的加、减、乘、除、乘方是数组的每个元素对该标量施加相应的加、减、乘、除、乘方设：a=[a1,a2,…,an],c=标量则：a+c=[a1+c,a2+c,…,an+c]
a.*c=[a1*c,a2*c,…,an*c]
a./c= [a1/c,a2/c,…,an/c](右除)
a.\c= [c/a1,c/a2,…,c/an] (左除)
a.^c= [a1^c,a2^c,…,an^c]
c.^a= [c^a1,c^a2,…,c^an]
程序,
a=[1 2 3 4]；
c=2；
a1=a+c
a2=a*c
a3=a./c
a4=a.\c
a5=a.^c
a6=c.^a
结果：
a1=3 4 5 6
a2=2 4 6 8
a3=0.5000 1.0000 1.5000 2.0000
a4=2.0000 1.0000 0.6667 0.5000
a5=1 4 9 16
a6=2 4 8 16
（2） 数组-数组运算
当两个数组有相同维数时，加、减、乘、除、幂运算可按元素对元素方式进行的，不同大小或维数的数组是不能进行运算的.
设：a=[a1,a2,…,an],b=[b1,b2,…,bn]
则：a+b= [a1+b1,a2+b2,…,an+bn]
a.*b= [a1*b1,a2*b2,…,an*bn]
a./b= [a1/b1,a2/b2,…,an/bn]
a.\b=[b1/a1,b2/a2,…,bn/an]
a.^b=[a1^b1,a2^b2,…,an^bn]
程序,
a=[2 2 2]；
b=[3 3 3]；
c1=a+b
c2=a.*b
c3=a./b
c4=a.\b
c5=a.^b
结果：
c1=5 5 5
c2=6 6 6
c3=0.6667 0.6667 0.6667
c4=1.5000 1.5000 1.5000
c5=8 8 8
§1.3 矩阵
§1.3.1 矩阵的建立
逗号或空格用于分隔某一行的元素，分号用于区分不同的行,除了分号，在输入矩阵时，按Enter键也表示开始一新行,输入矩阵时，严格要求所有行有相同的列,
例1.3 m=[1 2 3 4 ；5 6 7 8；9 10 11 12]
p=[1 1 1 1
2 2 2 2
3 3 3 3]
特殊矩阵的建立：
a=[ ] 产生一个空矩阵，当对一项操作无结果时，返回空矩阵，空矩阵的大小为零；
b=zeros(m，n) 产生一个m行、n列的零矩阵；
c=ones(m，n) 产生一个m行、n列的元素全为1的矩阵；
d=eye(m，n) 产生一个m行、n列的单位矩阵。
程序,
m=[1 2 3 4;5 6 7 8;9 10 11 12]
p=[1 1 1 1
2 2 2 2
3 3 3 3]
a=[]
b=zeros(2,3)
c=ones(2,3)
d=eye(2,3)
e=eye(3,3)
结果：
m =1 2 3 4
5 6 7 8
9 10 11 12
p =1 1 1 1
2 2 2 2
3 3 3 3
a =[]
b =0 0 0
0 0 0
c =1 1 1
1 1 1
d =1 0 0
0 1 0
e =1 0 0
0 1 0
0 0 1
§1.3.2 矩阵的运算
（1）标量-矩阵运算同标量-数组运算。
（2）矩阵-矩阵运算
[1] 元素对元素的运算，同数组-数组运算。
[2]矩阵运算：
矩阵加法：A+B
矩阵乘法：A*B
方阵的行列式：det（A）
方阵的逆：inv（A）
方阵的特征值与特征向量：[V，D]=eig[A]
程序,
a=[1 2 3
4 5 6]；
b=[1 2
1 2
1 2]；
c1=a+a
c2=a*b
c=[2 7 3;3 9 4;1 5 3]
c3=det(c)
c4=inv(c)
[v,d]=eig(c)
结果：
c1 =2 4 6
8 10 12
c2 = 6 12
15 30
c =2 7 3
3 9 4
1 5 3
c3 =-3
c4 =-2.3333 2.0000 -0.3333
1.6667 -1.0000 -0.3333
-2.0000 1.0000 1.0000
v =-0.5515 -0.7857 -0.2743
-0.7309 0.4412 -0.3391
-0.4020 -0.4337 0.8999
d =3.4635 0 0
0 -0.2747 0
0 0 0.8112
§1.4 关系与逻辑运算
§1.4.1 关系操作符
关系操作符
说明
<
小于
<=
小于或等于
>
大于
>=
大于或等于
= =
等于
~=
不等于
§1.4.2 逻辑运算符逻辑操作符
说明
＆
与
︱
或
~
非
§1.4.3 控制流
MATLAB提供三种决策或控制流结构：for循环、while循环、if-else-end结构,
这些结构经常包含大量的MATLAB命令，故经常出现在MATLAB程序中，而不是直接加在MATLAB提示符下.
a、for循环：
允许一组命令以固定的和预定的次数重复
for x=array
{commands}
end
在for和end语句之间的命令串{commands}按数组（array）中的每一列执行一次,在每一次迭代中，x被指定为数组的下一列，即在第n次循环中，x=array(：，n)
例1.4 对，求的值；
程序,
for n=1:10
x(n)=sin(n*pi/10);
end
x
结果：
x =0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000
b、While循环与for循环以固定次数求一组命令相反，while循环以不定的次数求一组语句的值.
while expression
{commands}
end
只要在表达式(expression)里的所有元素为真，就执行while和end语句之间的命令串{commands},
例1.5 设银行年利率为11.25%。将10000元钱存入银行，问多长时间会连本带利翻一番？
程序,
money=10000；years=0；
while money<20000
years=years+1；money=money*(1+11.25/100)；
end
years
money
结果：
years=7
money=2.1091e+004
c、If-Else-End结构
（1）有一个选择的一般形式是：
if expression
{commands}
end
如果在表达式(expression)里的所有元素为真，就执行if和end语句之间的命令串{commands}。
例 1.6编写程序 求函数 对任意的输入所得到的函数值。
先建立M文件fun1.m定义函数f（x），再在Matlab命令窗口输入fun1(2),fun1(-1)即可。
程序,
function f=fun1(x)
if x>1
f=x^2+1
end
if x<=1
f=2*x
end
输入：fun1(2),fun1(-1)
结果：5 -2
（2) 有三个或更多的选择的一般形式是：
if （expression1）
{commands1}
else if （expression2）
{commands2}
else if （expression3）
{commands3}
else if
……
else
{commands}
end
end
end
end
例1.7编写程序求函数 对任意的输入所得到的函数值。
先建立M文件fun2.m定义函数f（x），再在Matlab命令窗口输入fun2(2),fun2(0.5),fun2(-1)即可。
程序,
function f=fun2(x)
if x>1
f=x^2+1
else if x<=0
f=x^3
else
f=2*x
end
end
输入：
fun2(2),fun2(0.5),fun2(-1)
结果：
5 1 -1
§2 MATLAB作图
§2.1 二维图形
§2.1.1 曲线图
Matlab作图是通过描点、连线来实现的，故在画一个曲线图形之前，必须先取得该图形上的一系列的点的坐标（即横坐标和纵坐标），然后将该点集的坐标传给Matlab函数画图,命令为：
plot(X,Y,S)plot(X,Y)--画实线
plot(X,Y1,S1,X,Y2,S2,……,X,Yn,Sn)--将多条线画在一起X,Y是向量,
分别表示点集的横坐标和纵坐标，
plot(X,Y,S)描绘该点集所表示的曲线，其线型S确定如下：
y 黄色,点 - 连线
m 洋红 o 圈,短虚线
c 蓝绿色 x x-符号 -,长短线
r 红色 + 加号 -- 长虚线例2.1在[0,2*pi]用红线画sin(x),用蓝色圈画cos(x).
程序：
x=linspace(0,2*pi,30);
y=sin(x);
z=cos(x);
plot(x,y,'r',x,z,'bo')
结果：
§2.1.2 符号函数(显函数、隐函数和参数方程)画图
(1) ezplot
ezplot(‘f(x,y)’,[xmin,xmax,ymin,ymax])
表示在区间xmin<x<xmax和 ymin<y<ymax绘制隐函数f(x,y)=0的函数图
ezplot(‘x(t)’,’y(t)’,[tmin,tmax])
表示在区间tmin<t<tmax绘制参数方程x=x(t),y=y(t)的函数图例2.2 在[0,pi]上画y=sin(x)的图形.
程序,
ezplot(‘sin(x)’,[0,pi])
结果：

例2.3 在[0,2*pi]上画，星形图.
程序,
ezplot(‘cos(t)^3’,’sin(t)^3’,[0,2*pi])
结果：

例2.4在[-2，0.5]，[0，2]上画隐函数图.
程序,
ezplot('exp(x)+sin(x*y)',[-2,0.5,0,2])
结果：

§2.2 三维图形空间曲线
§2.2.1 一条曲线
plot3(x,y,z,s) x,y,z是n维向量，分别表示曲线上点集的横坐标、纵坐标、函数值;s指定颜色、线型等例2.5 在区间[0，10*pi]画出参数曲线x=sin(t),y=cos(t),z=t.
程序,
t=0:pi/50:10*pi;
plot3(sin(t),cos(t),t)
rotate3d %旋转结果：
§2.2.2 多条曲线
plot3(x,y,z) 其中x，y，z是都是m*n矩阵，其对应的每一列表示一条曲线.
例2.6 画多条曲线观察函数Z=(X+Y).^2,
程序：
x=-3:0.1:3;y=1:0.1:5;
[X,Y]=meshgrid(x,y);
Z=(X+Y).^2;
plot3(X,Y,Z)
结果：
（这里meshgrid(x,y)的作用是产生一个以向量x为行、向量y为列的矩阵）
§2.3 三维空间曲面
(1) surf(x,y,z)：画出数据点（x，y，z）表示的曲面。
例2.7 画函数Z=(X+Y).^2的图形,
程序：
x=-3:0.1:3;
y=1:0.1:5;
[X,Y]=meshgrid(x,y);
Z=(X+Y).^2;
surf(X,Y,Z)
shading flat %将当前图形变得平滑结果：
（2） Mesh(x,y,z)：画网格曲面。
例2.8 画出曲面Z=(X+Y).^2在不同视角的网格图程序：
x=-3:0.1:3; y=1:0.1:5;
[X,Y]=meshgrid(x,y);
Z=(X+Y).^2;
mesh(X,Y,Z)
结果：
(3) meshz(X,Y,Z) 在网格周围画一个curtain图(如,参考平面)。
例2.9 绘peaks的网格图程序：
[X,Y]=meshgrid(-3:.125:3);
Z=peaks(X,Y);
Meshz(X,Y,Z)
结果：
§2.4 图形处理
§2.4.1 图形上加格栅、图例和标注
（1）grid on,加格栅在当前图上；
grid off,删除格栅；
（2）hh = xlabel(string),在当前图形的x轴上加图例string；
hh = ylabel(string),在当前图形的y轴上加图例string；
hh = zlabel(string),在当前图形的z轴上加图例string；
hh = title(string),在当前图形的顶端上加图例string。
例2.10 在区间[0,2*pi]画sin(x)的图形，并加注图例“自变量X”、“函数Y”、“示意图”,并加格栅。
程序：
x=linspace(0,2*pi,30);y=sin(x);
plot(x,y)
xlabel('自变量X')
ylabel('函数Y')
title('示意图')
grid on
结果：
（3）hh = gtext(‘string’)运行命令gtext(‘string’)时，屏幕上出现当前图形，在图形上出现一个交叉的十字，该十字随鼠标的移动而移动，当按下鼠标左键时，该标注string放在当前十交叉的位置。
例2.11 在区间[0,2*pi]画sin(x)，并分别标注“sin(x)””cos(x)”.
程序：
x=linspace(0,2*pi,30);y=sin(x);z=cos(x);
plot(x,y,x,z)
gtext(‘sin(x)’);gtext(’cos(x)’)
结果：
0
§2.4.2 定制坐标
Axis([xmin xmax ymin ymax zmin zmax]) 定制图形坐标
Axis auto 将坐标轴返回到自动缺省值例2.12 在区间[0.005,0.01]显示sin(1/x)的图形。
程序：
x=linspace(0.0001,0.01,1000);y=sin(1./x);
plot(x,y)
axis([0.005 0.01 –1 1])
结果：
§2.4.3 图形保持
(1) hold on 保持当前图形,以便继续画图到当前图上。
hold off 释放当前图形窗口。
例2.13 将y=sin(x)、y=cos(x)分别用点和线画出在同一屏幕上。
程序：
x=linspace(0,2*pi,30);
y=sin(x);
z=cos(x)
plot(x,z,:)
hold on
plot(x,y)
结果：
(2) figure(h)
新建h窗口，激活图形使其可见，并把它置于其它图形之上例2.14 区间[0,2*pi]新建两个窗口分别画出y=sin(x)；z=cos(x)
程序：
x=linspace(0,2*pi,100);
y=sin(x);z=cos(x);
plot(x,y);
title('sin(x)');
pause
figure(2);
plot(x,z);
title('cos(x)');
结果：
§2.4.4 分割窗口
h=subplot(mrows,ncols,thisplot)
划分整个作图区域为mrows*ncols块（逐行对块访问）并激活第thisplot块，其后的作图语句将图形画在该块上。
例2.15将屏幕分割为四块，并分别画出y=sin(x)，z=cos(x)，a=sin(x)*cos(x),b=sin(x)/cos(x)。
程序：
x=linspace(0,2*pi,100);
y=sin(x); z=cos(x);
a=sin(x).*cos(x);b=sin(x)./(cos(x)+eps)
subplot(2,2,1);plot(x,y),title(‘sin(x)’)
subplot(2,2,2);plot(x,z),title(‘cos(x)’)
subplot(2,2,3);plot(x,a),title(‘sin(x)cos(x)’)
subplot(2,2,4);plot(x,b),title(‘sin(x)/cos(x)’)
结果：
§2.4.5 改变视角view
（1）view(a,b)：命令view(a,b)改变视角到(a,b),a是方位角,b为仰角。缺省视角为（-37.5，30）。
（2）view（[x，y，z]）：view用空间矢量表示的，三个量只关心它们的比例，与数值的大小无关，x轴view（[1，0，0]），y轴view（[0，1，0]），z轴view（[0，0，1]）。
例2.16 画出曲面Z=(X+Y).^2在不同视角的网格图.
程序：
x=-3:0.1:3; y=1:0.1:5;
[X,Y]=meshgrid(x,y);
Z=(X+Y).^2;
subplot(2,2,1); mesh(X,Y,Z)
subplot(2,2,2);mesh(X,Y,Z);view(50,-34)
subplot(2,2,3);mesh(X,Y,Z);view(-60,70)
subplot(2,2,4);mesh(X,Y,Z);view([0,1,1])
结果：
§2.4.6 动画制作
第一种方法是先保存多副不同的图片，然后连续回播；Moviein(),getframe,movie() 函数Moviein()产生一个帧矩阵来存放动画中的帧；函数getframe对当前的图象进行快照；函数movie()按顺序回放各帧。
例2.17将曲面peaks做成动画。
程序：
[x,y,z]=peaks(30);
surf(x,y,z)
axis([-3 3 -3 3 -10 10])
m=moviein(15);
for i=1:15
view(-37.5+24*(i-1),30)
m(:,i)=getframe;
end
movie(m)
结果：
说明：取所有图象中的两帧第二种方法是连续不段地擦除并重画屏幕对象，在重画时屏幕对象不断变化，实现动画效果，适于表现精度不高能够快速重画的场合。
程序：
s=0.05;
n=input('please input the number n=:')
x=rand(n,1)-0.5;y=rand(n,1)-0.5;
h=plot(x,y,'.');
axis([-2 2 -2 2])
set(h,'EraseMode','xor','MarkerSize',12)
i=1;pause
while i<n
drawnow
x=x+s*(-1)^i*rand(n,1);y=y+s*(-1)^i*rand(n,1);
set(h,'XData',x,'YData',y)
i=i+1;
end
输入：400
结果：<在动画始末时的两幅图象>
§2.5 特殊二、三维图形极坐标图：polar (theta,rho,s) 用角度theta（弧度表示）和极半径rho作极坐标图，用s指定线型。
例2.18 的极坐标图象程序：
theta=linspace(0,2*pi),
rho=sin(2*theta).*cos(2*theta);
polar(theta,rho,’g’)
title(‘Polar plot of sin(2*theta).*cos(2*theta)’);
结果：
第三篇 数学分支中的相关数学模型
§1 高等数学相关模型
§1.1 卫星轨道长度
人造地球卫星轨道可视为平面上的椭圆.我国第一颗人造地球卫星近地点距地球表面439km，远地点距地球表面2384km，地球半径为6371km,求该卫星的轨道长度.
分析：
卫星轨道椭圆的参数方程为分别是长、短半轴.根据计算参数方程的弧长公式，椭圆长度可表为如下积分
 （3-1-1）
称为椭圆积分，它无法用解析方法计算.
根据所给数据.
下面是用梯形公式和辛普森公式计算的程序：
function y=x5(t)
a=8755;b=6810;
y=sqrt(a^2*sin(t).^2+b^2*cos(t).^2);
在MATLAB工作区键入以下程序：（为区别两段程序，以下用空行隔开）
t=0:pi/10：pi/2
y1=x5(t);
L1=4*trapz(t,y1)
L2=4*quad(‘x5’,0,pi/2,le-6)
用长格式输出结果为：
L1=4.908996526785276e+004
L2=4.908996531830460e+004
可以看出，梯形公式（仅将区间5等分）给出了很好的结果.轨道长度为4.90910千米.
§1.2 射击命中概率
炮弹射击的目标为一正椭圆形区域，当瞄准目标的中心发射时，在众多因素的影响下，弹着点与目标中心有随机偏差.可以合理的假设弹着点围绕中心成二维正态分布，且偏差在方向和方向相互独立.若椭圆形区域在方向半轴长120m，方向半轴长80，设弹着点偏差的均方差在方向和方向均为100m,求炮弹落在椭圆形区域内的概率.
分析：
按题目所述，弹着点与目标中心的偏差服从二维正态分布，且在方向和方向相互独立.设目标中心为则弹着点的概率密度函数是
 （3-1-2）
其中(m),炮弹命中椭圆形区域的概率为二重积分
 （3-1-3）
其中(m).它无法用解析方法求解，我们用蒙特卡罗方法计算.
作变换且以100(m)为1单位，即，于是

 （3-1-4）
问题已化为蒙特卡罗方法计算二维积分的形式，可以编制以下程序计算：
a=1.2;b=0.8;m=0;z=0;n=100000;
for i=1:n
x=rand(1,2);
y=0;
if x(1)^2+x(2)^2<=1
y=exp(-0.5*(a^2*x(1)^2+b^2*x(2)^2));
z=z+y
m=m+1
end
end
p=4*a*b*z/2/pi/n,m
这里次，结果为命中概率.实际上，结果具有随机性，不过概率在0.37~0.38是可信的.
课堂练习：冰淇淋的下部为锥体，上部为半球，设它由锥面和球面围成，
（1）基于适当假设建立求冰淇淋体积的数学模型；
（2）对于你所建立的模型，能否求其精确解？如果能，求其精确解，如果不能，试基于Monte-carlo（蒙特卡罗）方法写出在matlab中可实现模型求解的程序；
班级,姓名,学号,成绩,
§1.3 人口增长率
已知20世纪美国人口统计数据如下，试计算表3-1-1中这些年份的人口增长率.
表3-1-1 20世纪美国人口统计数据年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990
人口（） 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4
又已知某地区20世纪70年代的人口增长率如表3-1-2，且1970年人口为210（百万），试估计1980年的人口.
表3-1-2 某地区20世纪70年代人口增长率数据年份 1970 1972 1974 1976 1978 1980
年增长率（%） 0.87 0.85 0.89 0.91 0.95 1.10
分析：
（1）若记时刻的人口为，则人口（相对）增长率为表示每年人口增长的比例.对于题目给出的人口数据（表3-1-2），记1900年为
年依次为相应地，人口记为，年增长率为，用数值微分的三点公式计算（将每10年的增长率变为每年的增长率）


结果为,
年 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990
 2.20 1.66 1.46 1.02 1.04 1.58 1.49 1.16 1.05 1.04
可以看出，20世纪美国人口增长率总的来说在下降，但是有起伏，30年代和二次大战时期人口增长率显著下降，战后又迅速上升，80年代和90年代稳定在1%多一点.
（2）仍用上面的记号，人口增长满足微分方程它在初始条件下的解为

因为在题目中增长率以离散数据给出（表3-1-2），所以上式要用数值积分计算。用梯形公式计算的结果是，1980年该地区人口为230.2（百万）。
附,数值微分的三点公式：

思考题,国土面积问题为了算出瑞士的国土面积，首先对瑞士地图作如下测量：以由西向东方向为轴，由南到北方向为轴，选择方便的原点，并将从最西边界点到最东边界点在轴上的区间适当地划分为若干段，在每个分点的方向测出南边界点和北边界点的坐标和，这样就得到了表3-1-3中的数据（单位mm）.
根据地图的比例我们知道18mm相当于40km，试由测量数据计算瑞士国土的近似面积，与它的精确值41288km比较.
表3-1-3 瑞士地图测量数据

7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80.5 91.0

44 45 47 50 50 3 8 30 30 34 36 34 41 45 46

44 59 70 72 93 100 110 110 110 117 118 116 118

96 101 104 106.5 111.5 118 123.5 136.5 142 146 150 157 158

43 37 33 28 32 65 55 54 52 50 66 66 68

121 124 121 121 121 122 116 83 81 82 86 85 68
§2 线性代数相关模型
§2.1 投入产出综合平衡分析
国民经济各个部门之间存在着相互依存的关系，每个部门在运转中将其他部门的产品或半成品经过加工（称为投入）变为自己的产品（称为产出），如何根据各部门间的投入—产出关系，确定各部门的产出水平，以满足社会的需求，是投入产出综合平衡模型研究的课题.试讨论如下的简化问题.
设国民经济仅由农业、制造业、和服务业三个部门构成，已知某年它们之间的投入产出关系、外部需求、初始投入等如表3-2-1所示（数字表示产值，单位为亿元）.
表3-2-1 国民经济各个部门间的关系产出投入
农业
制造业
服务业
外部需求
总产出
农业
15
20
30
35
100
制造业
30
10
45
115
200
服务业
20
60
/
70
150
外部需求
35
110
75
总投入
100
200
150
表中第一行数字表示，农业总产出为100亿元，其中15亿元农产品用于农业生产本身（如提供种子），20亿元用于制造业（如提供木材、毛皮），30亿元用于服务业，剩下35亿元农产品用来满足外部需求（包括消费、积累、出口等）.可以类似的解释第二、三、行数字.第一列数字中，15亿元如前所述，30亿元是制造业对农业的投入（如提供农具），20亿元是服务业对农业的投入，35亿元的初始投入包括工资、税收、进口等，总投入100亿元与总产出相等.
假定每个部门的产出与投入是成正比的，由表3-2-1能够确定这三个部门的投入产出表，如表3-2-2。
表3-2-2 投入产出表
产出投入
农业
制造业
服务业
农业
0.15
0.10
0.20
制造业
0.30
0.05
0.30
服务业
0.20
0.30
0
表中第一行、第二列的数字0.10表示，生产1个单位产值的制造业产品需投入0.10个单位产值的农产品，这是由表1中20亿元农产品投入制造业，可以产出200亿元制造业总产值而来的（20200=0.1）.同样，第三行、第一列的数字0.20表示，生产1个单位产值的农产品需要0.20个单位的服务业产值，因为表3-2-1中20亿元的服务业产值投入农业，得到100亿元的农业总产值.表3-2-2的数字称为投入系数或消耗系数，如果技术水平没有变化，可以假设系数是常数.
（1）设有个部门，已知投入系数，给定外部需求，建立求解各部门总产出的模型.
（2）系数如表3-2-2所给，如果今年对农业、制造业、和服务业的外部需求分别为50，150，100亿元，问这三个部门的总产出分别应为多少？
（3）如果三个部门的外部需求分别增加1个单位，他们的总产出应分别增加多少？
（4）如果对于任意给定的、非负的外部需求，都能得到非负的总产出，模型就称为可行的.问为使模型可行，投入系数应满足什么条件？
分析：
（1） 投入产出综合平衡分析
① 若共有个部门，记一定时期内第个部门的总产出为，其中对第个部门的投入为，满足的外部需求为，则
  （3-2-1）
表3-2-1的每一行既满足（3-2-1）式.记第个部门的单位产出需要第个部门的投入为，在每个部门的产出与投入成正比的假定下，有
  （3-2-2）表5中的投入系数即为,（3-2-2）式代入（3-2-1）式得
  （3-2-3）
记投入系数矩阵，产出向量，需求向量，则（3-2-3）式写作
 （3-2-4）
或
 （3-2-5）若可逆，则
 （3-2-6）
当投入系数和外部需求给定后，即可算出各部门的总产出.
将表3-2-2中的投入系数及对三个部门的外部需求输入MATLAB，用以下程序：
a=[0.15 0.1 0.2;0.3 0.05 0.3;0.2 0.3 0];
d=[50 150 100];
b=eye(3)-a;
x=b\d
c=inv(b)
得到三个部门的总产出分别应为139.2801，267.6056，208.1377（亿元）
③（3-2-6）式表明总产出对外部需求是线性的，所以当增加1个单位时，的增量由决定.在上面程序中已经算出为
1.3459 0.2504 0.3443
0.5634 1.2676 0.4930
0.4382 0.4304 1.2167
其第一列数字表明，当对农业的需求增加1个单位时，农业、制造业、和服务业的总产出应分别增加1.3459，0.5634，0.4382单位.其余类似，这些数字称部门关联系数。
④ 为了使模型可行，即对任意的外部需求，由（3-2-6）式能够得到总产出，显然只需（指每个元素非负，下同）
因为（由定义（3-2-2）式），

所以只要就有.由矩阵范数定义可知与等价，且，故只要（这里取便于应用的1—范数），即
  （3-2-7）
模型就是可行的.
如果投入系数是根据实际数据算出的（如由表3-2-1得到表3-2-2），由（3-2-2）式可知（3-2-7）式等价于
  （3-2-8）可以看出，只要初始投入非负，（3-2-8）式是自然成立的。
课堂练习：现有一个木工、一个电工和一个油漆工，三人相互彼此装修他们自己的房子，在装修之前，他们达成了如下协议：（1）每人总共工作10天（包括给自己家干活在内）；（2）每人的日工资根据一般的市价在60~80元之间；（3）每人的日工资数应使得每人的总收入与总支出相等。表1是他们协商后制定出的工作天数的分配方案,
表1
工种天数
木工
电工
油漆工
在木工家的工作天数
2
1
6
在电工家的工作天数
4
5
1
在油漆工家的工作天数
4
4
3
（1）、建立相应的数学模型计算出他们每人应得的工资；
（2）、对于一般情形由个人构成的相互装修的系统，给出他们应得工资的模型，并编写出计算机上可实现的程序；
班级,姓名,学号,成绩,
§2.2 输电网络
一种大型输电网络可简化为图3-2-1所示电路，其中表示负载电阻表示线路内阻，设电源电压为.
（1）列出求各负载上电流的方程；
  
     

图3-2-1 输电网络简化电路
（2）设在
的情况下求及总电流；
（3）讨论的情况.
分析：
（1）记上的电流为根据电路中电流、电压的关系可以列出
 （3-2-9）
和
 （3-2-10）
消去得到
 （3-2-11）
记
 
则方程（3-2-11）表为
 （3-2-12）
（3-2-11）或（3-2-12）即为求解电流的方程.
（2）将已知条件输入，编写以下程序，并注意到总电流为各负载上电流之和.
r=1;R=6;v=18;n=10;
b1=sparse(1,1,v,n,1);
b=full(b1);
a1=triu(r*ones(n,n));
a2=diag(R*ones(1,n));
a3=-tril(R*ones(n,n),-1)+ tril(R*ones(n,n),-2);
a=a1+a2+a3;
I=a\b;
I0=sum(I)
计算结果如表3-2-3所示（作为比较，表中还列出了的结果）：
表3-2-3 输电网络电流计算结果

0
1
2
3
4
5
6
7
8
9
10

5.9970
6.0000
2.0005
2.0000
1.3344
1.3333
0.8907
0.8889
0.5955
0.5926
0.3995
0.3951
0.2702
0.2634
0.1858
0.1756
0.1324
0.1171
0.1011
0.0780
0.0867
0.0520

11
12
13
14
15
16
17
18
19
20

0.0347
0.0231
0.0154
0.0103
0.0069
0.0047
0.0032
0.0023
0.0018
0.0015
可以看出，从到，总电流几乎不变，变化也很小，
并且差不多是的倍.如果增加到50，100，可以得到类似的结论.为证实这一结论需要研究的情况.
（3）上面的计算看出，增大时总电流几乎不变，说明总电阻增加很少.为求出总电阻考察图3-2-2所示的第段电路（从右向左数）和第段电路的等效电阻和，有
 （3-2-13）
以代入可计算，当时记为，满足

即

由此解得
  （3-2-14）

  
图3-2-2 输电网络的等效电阻
于是总电流

而
 

不难推出
 
证明了前面的计算结果.
习题：种群的繁殖与稳定收获种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应维持不变.种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性.
种群年龄记作当年年龄的种群数量记作，繁殖率记作（每个雌性个体一年繁殖的数量），自然存活率记作为一年的死亡率），收获量记作，则来年年龄的种群数量应为

（1）若已知，给定收获量，建立求各年龄的稳定种群数量的模型（用矩阵、向量表示）。
（2）设如要求为500，400，200，100，100，求。
（3）使均为500，如何达到?
§3 概率统计相关模型
§3.1合金强度与碳含量
合金的强度与其中的碳含量有比较密切的关系，今从生产中收集了一批数据如表3-3-1，试先拟合一个函数，再用回归分析对它进行检验.
表3-3-1 合金的强度（kg/mm）与其中的碳含量

0.10
0.11
0.12
0.13
0.14
0.15
0.16
0.17
0.18
0.20
0.21
0.23

42.0
41.5
45.0
45.5
45.0
47.5
49.0
55.0
50.0
55.0
55.5
60.5
将的12个点画在图3-3-1上，可知与大致上为线性关系，用polyfit(x,y,1)拟合得到，图3-3-1也给出了拟合的直线.
程序：
x=0.1:0.01:0.23;
x=[x(1:9),x(11:12),x(14)];
y=[42,41.5,45,45.5,45,47.5,49,55,50,55,55.5,60.5];
pp=polyfit(x,y,1);
xx=0.08:0.01:0.25;
yy=polyval(pp,xx);
plot(x,y,'r*',xx,yy)
结果：

图3-3-1 合金强度与碳含量拟合图这个实验要对上述结果的可信度进行检验，判断对的影响是否显著，检查数据中有无异常点（0.17，55.0）是可疑的异常点，由的取值对做出预测等.
分析：先介绍MATLAB统计工具箱
MATLAB统计工具箱用命令regress实现多元线性回归，用法是：
b=regress(Y,X) Y,X为按列排列的数据，b为回归系数估计值
[b,bint,r,rint,stats]=regress(Y,X,alpha) Y,X同上，alpha为显著性水平（缺省时设定为0.05），b,bint为回归系数估计值和他们的置信区间，3个数值，第1个是，第2个是，第3个是与对应的概率时拒绝，回归模型成立,残差及其置信区间可以用rcoplot(r,rint)画图.
注：其中称为样本方差；，称为回归平方和，称为残差平方和；，越大，与相关关系越密切，
合金强度与碳含量问题回归模型为
 （3-3-1）
输入表3-3-1数据，用regess和rcoplot编程：
x1=0.1:0.01:0.18;
x=[x1 0.2 0.21 0.23]’;
y=[42 41.5 45 45.5 45 47.5 49 55 50 55 55.5 60.5]’;
x=[ones(12,1) x];
[b,bint,r,rint,stats]=regress(y,x);
b,bint,stats,rcoplot(r,rint)
得到
b=27.0269 140.6194
bint=22.3226 31.7313
111.7842 169.4546
stats=0.9219 118.0670 0.0000
即 
的置信区间是[22.3226，31.7313]，的置信区间是[111.7842，169.4546]；
可知模型（3-3-1）成立.
观察图3-3-2所示的残差分布，除第8个数据外其余残差的置信区间均包含零点，第8个点应视为异常点，将其剔除后重新计算，可得：
b=26.8968 139.9043
bint=24.1330 29.6606
122.7939 157.0148
stats=0.9744 342.1259 0.0000
我们看到，，变化不大，但其置信区间明显缩小，和都明显变大，应该用修改后的这个结果.

图3-3-2 合金强度与碳含量的残差图
§3.2 年龄与运动能力
如果从数据的散点图上发现与成较明显的二次（或高次）函数关系，或者用线性模型（3-3-1）的效果不太好，就可以选用多项式回归，见下面的例子：
将17至29岁的运动员每两岁一组分为7组，每组两人测量其旋转定向能力，以考察年龄对这种运动能力的影响.现得到一组数据如表3-3-2，试建立二者之间的关系.
表3-3-2 年龄与旋转定向能力的测量数据年龄
17
19
21
23
25
27
29
第一人
20.48
25.13
26.15
30.0
26.1
20.3
19.35
第二人
24.35
28.11
26.3
31.4
26.92
25.7
21.3
数据的散点图如图3-3-3，明显地呈现两端低中间高的形状，所以应拟合一条二次曲线.虽然用ployfit能完成这一任务，但仍需作进一步的统计分析.我们还将介绍其它更合适的命令.
程序：x=17:2:29;
y1=[20.48,25.13,26.15,30.0,26.1,20.3,19.35];
y2=[24.35,28.11,26.3,31.4,26.92,25.7,21.3];
plot(x,y1,'+',x,y2,'+')
axis([15 30 15 35])
结果：

图3-3-3 年龄与旋转定向能力
§3.2.1一元多项式回归
一元多项式回归可用命令polyfit实现，我们通过年龄与运动能力问题说明其用法.选用二次模型
 （3-3-2）
输入表3-3-2的数据并编程：
x1=17:2:29;x=[x1,x1];
y=[20.48 25.13 26.15 30.0 26.1 20.3 19.35…24.35 28.11 26.3 31.4 26.92 25.7 21.3];
[p,s]=polyfit(x,y,2);
p
得到
p=-0.2003 8.9782 -72.2150
即
a1=-0.2003 a2=8.9782 a3=-72.2150.
上面的S是一个数据结构，用于其它函数的计算，如
[Y,delta]=polyconf(p,x,s);Y
得到的拟合值（只列出前7个，后7个与此相同）
Y=22.5243 26.0582 27.9896 28.3186 27.0450 24.1689 19.6904
作与的（连续）曲线，拟合的效果从图3-3-4可以直观地看出.

图3-3-4 年龄与运动能力的拟合曲线
 和可作为衡量拟合的优劣的指标，作如下计算：
y1=mean(y);
resquare=sum((Y-y1).^2)./sum((y-y1).^2),
s=sqrt(sum((y-Y).^2)./12),
得到
resquare=0.6980
s=2.0831
这是一个尚能接受的结果

图3-3-5 年龄与运动能力拟合的交互式界面用polytool(x,y,2)可以得到一个如图3-3-5的交互式画面，在画面中绿色曲线为拟合曲线，它两侧的红线是的置信区间.你可以用鼠标移动图中的十字线来改变图下方的值，也可以在窗口内输入，左边就给出的预测值及其置信区间.通过左下方的Export下拉式菜单，可以输出回归系数等，你不妨试试看.这个命令的用法与下面将介绍的rstool相似.
§3.3 商品销售量与价格
某厂生产的一种电器的销售量与竞争对手的价格和本厂的价格有关.表3-3-3是该商品在10个城市的销售记录，试根据这些数据建立与和的关系式，对得到的模型和系数进行检验.若某市本厂产品售价160（元），竞争对手售价170（元），预测商品在该市的销售量.
表3-3-3 商品销售量与价格和
（元）
120
140
190
130
155
175
125
145
180
150
（元）
100
110
90
150
210
150
250
270
300
250
（个）
102
100
120
77
46
93
26
69
65
85
将各10个点画在图3-3-6上，可以看出与有较明显的线性关系，而与之间的关系则难以确定，我们将作几种尝试，用统计分析决定优劣.

图3-3-6 商品销售量与价格和
分析,商品销售量与价格问题回归模型为
 （3-3-3）
输入表3-3-3数据并编程：
x1=[120 140 190 130 155 175 125 145 180 150];
x2=[100 110 90 150 210 150 250 270 300 250];
y=[102 100 120 77 46 93 26 69 65 85];
x=[ones (10,1)x1’x2’];
[b,bint,r,rint,stats]=regress(y,x);
得到
b=66.5176 0.4139 –0.2698
bint=-32.5060 165.5411
-0.2018 1.0296
-0.4611 -0.0785
stats=0.6527 6.5786 0.0247
可以看出结果不是太好：，取时回归模型（3-3-3）可用，
但取则模型不能用；，较小；，的置信区间包含了零点.下面将试图用的二次函数改进它.
统计工具箱提供了一个作多元二项式回归的命令rstool，它也产生一个交互式界面，并输出有关信息，用法是
rstool(x,y,model,alpha) 输入数据x,y分别为矩阵和维向量，alpha为显著性水平（缺省时设定为0.05），model由下列4个模型中选择1个（用字符串输入，缺省时设定为线性模型）：
linear（线性）：
purequadratic（纯二次）：
interaction（交叉）：
quadratic（完全二次）：
我们用它再做一遍商品销售量与价格问题，选择纯二次模型，即
 （3-3-4）
程序：
x1=[120,140,190,130,155,175,125,145,180,150];
x2=[100,110,90,150,210,150,250,270,300,250];
y=[102,100,120,77,46,93,26,69,65,85];
x=[x1’,x2’];
rstool(x,y,’ purequadratic’)
结果：

图3-3-7 商品销售量与价格拟合的交互式界面
得到一个如图3-3-7的交互式界面，给出两幅图形，左边是（=151）固定时的曲线及其置信区间，右边是（=188）固定时的曲线及其置信区间.用鼠标移动图中的十字线，或在图下方窗口内输入，可改变，.图左边给出的预测值及其置信区间，就用这种界面可以回答商品销售量与价格问题中提出的“若某市本厂产品售价160（元），竞争对手售价170（元），预测该市的销售量”问题.
图的左下方有两个下拉式菜单，一个菜单Export用以向MATLAB工作区传送数据，包括beta（回归系数），rmse（剩余标准差），residuals（残差）,模型（3-3-4）的回归系数和剩余标准差为
beta=-312.5871 7.2701 -1.7337 -0.0228 0.0037
rmse=16.6436
另一个菜单model用以在上述4个模型中选择，你可以比较一下它们的剩余标准差，会发现以模型（3-3-4）的rmse=16.6436最小.
附,逐步回归实际问题中影响因变量的因素可能很多，我们希望从中挑选出影响显著的自变量来建立回归模型，这就涉及到变量选择的问题，逐步回归是一种从众多变量中有效地选择重要变量的方法.以下只讨论线性回归模型式的情况.
变量选择的标准，简单地说就是所有对因变量影响显著的变量都应选入模型，而影响不显著的变量都不应选入模型，从便于应用的角度应使模型中变量个数尽可能减少.
若候选的自变量集合为，从中选出一个子集，设中有个自变量，由和因变量构造的回归模型的误差平方和，则模型的剩余标准差的平方 为数据样本容量.所选子集应使尽量小，通常回归模型中包含的自变量越多，误差平方和为越小，但若模型中包含有对影响很小的变量，那么不会由于包含这些变量在内而减少多少，却因的增加可能使最小作为衡量变量选择的一个数量标准。
逐步回归是实现变量选择的一种方法，基本思路为：先确定一初始子集，然后每次从子集外影响显著的变量中引入一个对影响最大的，再对原来子集中的变量进行检验，从变得不显著的变量中删除一个影响最小的，直到不能引入和剔除为止.使用逐步回归有两点值得注意，一是要适当地选定引入变量的显著性水平和删除变量的显著性水平，显然，越大，引入的变量越多；越大，剔除的变量越少.二是由于各个变量之间的相关性，一个新的变量引入后，会使原来认为显著的某个变量变得不显著，从而被剔除，所以在最初选择变量时应尽量选择相互独立性强的那些.
在MATLAB统计工具箱中用作逐步回归的是命令stepwise，它提供了一个交互式画面，通过这个工具你可以自由地选择变量，进行统计分析，其通常用法是：
stepwise(x,y,inmodel,alpha) x是自变量数据，y是因变量数据，分别为和矩阵，inmodel是矩阵x的列数的指标，给出初始模型中包括的子集（缺省时设定为全部自变量），alpha为显著性水平.
Stepwise命令产生三个图形窗口：Stepwise Regression Plot,Stepwise Regression Diagnostics Table,Stepwise History Plot.所有这些图形界面都有热点，即当鼠标移到图形的某个区域时，鼠标的指针会变成一个小圆，点击后会产生交互的作用.
在Stepwise Regression Plot窗口，显示回归系数及其置信区间，绿色表明在模型中的变量，红色表明从模型中移去的变量，两边有虚线或实线，虚线表示该变量的拟合系数与零无显著差异.在这个窗口中还有Scale Inputs 和Export按钮.
按下Scale Inputs表明对于输入数据的每列进行正态化处理，使其标准差为1.点击Export产生一个菜单，标明了要传送给MATLAB工作区的参数，他们给出了统计计算的一些结果.
Stepwise Regression Diagnostics Table窗口中列出了一个统计表，包括回归系数及其置信区间，模型的统计量（RMSE,R-square,F,p等，其含义与regress,rstool相同）.你可以通过这些统计量的变化来确定模型.
Stepwise History窗口显示RMSE的值及其置信区间，和每次选择不同模型后的结果.点击相关的区域就可以把对应的Stepwise Regression Plot窗口和Stepwise Regression Diagnostics Table窗口调出.下面通过一个例子说明stepwise的用法.
例 水泥固定时放出的热量与水泥中4种化学成分有关，今测得一组数据如下，使用逐步回归来确定一个线性模型。
序号





1
7
26
6
60
78.5
2
1
29
15
52
74.3
3
11
56
8
20
104.3
4
11
31
8
47
87.6
5
7
52
6
33
95.9
6
11
55
9
22
109.2
7
3
71
17
6
102.7
8
1
31
22
44
72.5
9
2
54
18
22
93.1
10
21
47
4
26
115.9
11
1
40
23
34
83.3
12
11
66
9
12
113.3
13
10
68
8
12
109.4
将数据按表中格式输入，x=[x1 x2 x3 x4],用
stepwise(x,y)
得到一个Stepwise Regression Diagnostics Table如下：
Confidence Intervals
Column Parameter Lower Upper
1 1.551 -0.8319 3.934
2 0.5102 -1.806 2.826
3 0.1019 -2.313 2.517
4 -0.1441 -2.413 2.125
RMSE R-square F P
2.446 0.9824 111.5 4.756e-007
可以看出，不显著，移去这两个变量后的统计结果如下：
Confidence Intervals
Column Parameter Lower Upper
1 1.468 1.1 1.836
2 0.6623 0.5232 0.8013
3 0.25 -0.3235 0.8236
4 -0.2365 -0.7746 0.3015
RMSE R-square F P
2.406 0.9787 229.5 4.407e-009
这个表中的两行用红色显示，表明它们已移去.
从新的统计结果可以看出，虽然剩余标准差(RMSE)没有太大的变化，但是统计量的值明显增大，因此新的回归模型更好一些，使用前面的回归分析方法可以求出最终的模型为

第四篇 典型案例分析
§1 投篮的出手角度
§1．1 问题的提出在激烈的篮球比赛中,提高投篮命中率对于获胜无疑起着决定作用,而出手角度和出手速度是决定投篮能否命中的两个关键因素。这里讨论比赛中最简单、但对于胜负也常常是很重要的一种投篮方式—罚球,并且球出手后不考虑自身的旋转,不考虑碰蓝板或篮筐。
图4-1-1为过罚点球和篮框中心、且垂直于地面的平面示意图,按照标准尺寸和点的水平距离m,点的高度m,篮球直径cm,篮球直径cm。不妨假定篮球运动员的出手高度为1.8~2.1m,出手速度为8.0~9.0m/s.试建立数学模型研究以下问题：
（1）先不考虑篮球和篮框的大小,讨论球心命中框心的条件。对不同的出手高度和出手速度,确定出手角度和篮框的入射角度;
（2）考虑篮球和篮框的大小,讨论球心命中框心且球入框的条件。检查上面得到的出手角度和入射角度是否符合这个条件；
（3）为了使球入框,球心不一定要命中框心,可以偏前或偏后(球心命中图6中的或点)。讨论保证球入框的条件下,出手角度允许的最大偏差,和出手速度允许的最大偏差；
（3）考虑空气阻力的影响。由于投篮基本上是水平方向且速度不大的室内运动,可以只计水平方向的阻力,设阻力与速度成正比,比例系数不超过0.05(1/s)。

 
P   
H
h
L
图4-1-1 从罚球点投篮示意图
§1．2 问题的分析
（1）不考虑篮球和篮框大小的简单情况,相当于将球视为质点(球心)的斜抛运动。将坐标圆点定在球心,列出(水平)方向和(竖直)方向的运动方程,就可以得到球心的运动轨迹,于是球心命中框心的条件可以表示为出手角度与出手速度、出手高度之间的关系,以及球框的入射角度与出手角度的关系,由此可对不同的出手速度和出手高度,计算出手角度和入射角度。
（2）考虑篮球和篮框的大小时.如图4-1-2,篮球直径为,篮框直径为.显然,即使球心命中框心，若入射角太小，球会碰到框的近侧，不能入框。由图4-1-2不难得出应满足的、球心命中框心且球入框的条件。前面计算结果中不满足这个条件的，当然应该去掉。
（3）球入框时球心可以偏离框心，偏前（图4-1-1的点）的最大距离为图4-1-3中的，可以从入射角算出。根据和球心轨迹中与的关系，能够得到出手角度允许的最大偏差。出手速度允许的最大偏差可以类似的处理。
（4）考虑水平方向的空气阻力时，应该用微分方程求解球心的运动轨迹，由于阻力很小，可作适当简化。然后与前面类似的作各种计算。
d
 O  d O
B A Δx B
D
D
图4-1-2 篮球入筐 图4-1-3 球心偏前
§1．3 基本模型
（1）不考虑篮球和篮框的大小，不考虑空气阻力的影响，以未出手时的球心为坐标原点，轴为水平方向，轴为竖直方向，篮球在时以出手速度和出手角度投出，可视为质点（球心）的斜抛运动，其运动方程是我们熟知的

 （4-1-1）
其中是重力加速度。由此可得球心运动轨迹为如下抛物线
 （4-1-2）
以代如（4-1-2）式，就得到球心命中框心的条件
 （4-1-3）
可以看出，给定出手速度和出手高度，有两个出手角度满足这个条件。而（4-1-3）式有解的前提为
 （4-1-4）
可对解得
 （4-1-5）
于是对于一定的出手高度，使（4-1-5）式等号成立的为最小出手速度，它是的减函数。由（4-1-3）式计算出的两个出手角度记作，，且设，可以看出，是和的增函数。
思考 是的减函数，是和的增函数，能做出实际解释吗？与和的关系又怎样呢？
（2）球入篮框时的入射角度可从下式得到
 （4-1-6）
这里的导数由（4-1-2）式计算代入后得
 （4-1-7）
思考 ，会小于零吗？会小于零吗？如果小于零，作何解释。
（3）考虑篮球和篮框的大小，如图4-1-2，若入射角度太小，则球无法入框。由图不难看出，球心命中框心且球入框的条件为
 （4-1-8）
将cm，cm代入得.
§1．4 出手角度和出手速度最大偏差估计由图4-1-3看出，球入框时球心可以偏前（偏后与偏前一样）的最大距离为
 （4-1-9）
为了得到出手角度允许的最大偏差，可以在（4-1-3）式中以代替重新计算，但是由于中包含，从而也包含，所以这种方法不能解析地求出。
如果从（4-1-2）式出发并将代入，可得
 （4-1-10）
对求导并令，就有
 （4-1-11）
用近似代替左边的导数，即可得到出手角度的偏差与的如下关系
 （4-1-12）
由和已经得到的也容易计算相对偏差。
类似地，（4-1-10）式对求导并令，可得出手速度允许的最大偏差
 （4-1-13）
由（4-1-12）、（4-1-13）式的相对偏差为
 （4-1-14）
§1．5 空气阻力的影响按照篮球运动的特点可以只考虑水平方向的阻力，且阻力与速度成正比，比例系数为。这时水平方向的运动由微分方程
   （4-1-15）
描述，其解为
  （4-1-16）
因为阻力不大（不超过0.05秒），时间也很小（约1秒），所以将（4-1-16）式中的作泰勒展开后忽略二阶以上项得到（不考虑竖直方向的阻力，故仍与（4-1-1）式相同）

 （4-1-17）
在不考虑篮球和篮框大小时，球心命中框心的条件由方程组

 （4-1-18）
确定。
§1．6 算法实现和计算结果
（1）对不同出手高度的最小出手速度和相应的出手角度使（4-1-5）式等号成立的为最小出手速度，在这个速度下由（4-1-3）式可得相应的出手角度为

取出手高度（m），计算结果见表4-1-1。
表4-1-1 对不同出手高度的最小出手速度和相应的出手角度
h(m)
(m/s)
(度)
1.8
7.6789
52.6012
1.9
7.5985
52.0181
2.0
7.5186
51.4290
2.1
7.4392
50.8344
（2）对不同出手速度和出手高度的出手角度和入射角度对出手速度（m/s）和出手高度（m），由（4-1-3）式计算出手角度由（4-1-7）式计算入射角度，结果见表4-1-2。
表4-1-2 对不同出手速度和出手高度的出手角度和入射角度
v（m/s）
h (m)
(度) (度) (度) (度)
8.0
1.8
1.9
2.0
2.1
62.4099 42.7925 53.8763 20.9213
63.1147 40.9188 55.8206 201431
63.7281 39.1300 57.4941 19.6478
64.2670 37.4017 5809615 19.3698
8.5
1.8
1.9
2.0
2.1
67.6975 37.5049 62.1726 12.6250
68.0288 36.0075 63.1884 12.7753
68.3367 34.5214 64.1179 13.0240
68.6244 33.0444 64.9279 13.3583
9.0
1.8
1.9
2.0
2.1
71.0697 34.1327 67.1426 7.6550
71.2749 3207614 67.7974 8.1663
71.4700 31.3881 68.4098 8.7321
71.6561 30.0127 68.9840 9.3472
可以看出，均小于，不满足（4-1-8）式的条件，所以在考虑篮球和篮框大小的实际情况下，出手角度只能是。
（3）出手角度和出手速度的最大偏差利用（4-1-12）式和上面的，计算出手角度最大偏差和，再用（4-1-13），（4-1-14）式计算出手速度的最大偏差和，只将=1.8，2.0（m）的结果列入表4-1-3。
表4-1-3出手角度和出手速度最大偏差
h(m)
(度) v(m/s)
 
 
1.8
62.4099 8.0
67.6975 8.5
71.0697 9.0
-0.7562 0.0528
-0.5603 0.0694
-0.4570 0.0803
1.2261 0.6597
0.8276 0.8167
0.6431 0.8925
2.0
63.7281 8.0
68.3367 8.5
71.4700 9.0
-0.7100 0.0601
-0.5411 0.0734
-0.4463 0.0832
1.1140 0.7511
0.7918 0.8640
0.6244 0.9243
（4）空气阻力的影响设空气阻力系数=0.05（1/s），对出手速度=8.0~9.0（m/s）和出手高度=1.8~2.1（m），由（4-1-18）式计算出手角度。（4-1-18）是非线性方程组，可用MATLAB优化工具包中的fsolve求解：
function f=qiu（x,v,g,k,H,h,L,）
f(1)=v*x(1)*x(2) -k*v* x(1)*x(2) ^2/2-L;%x(1)=cos,x(2)=t
f(2)=v*sqrt(1-x(1)^2)*x(2) -g*x(2) ^2/2-(H-h);
g=9.8;k=0.05;H=3.05;l=4.6;
x0=[0.4 1];%根据不考虑阻力的结果选取x的初值
for i=1:3
v(i)=8+(i-1)*0.5;
for j=1:4
h(j)=1.8+(j-1)*0.1;
x=fsolve(‘qiu’,x0,[],[],v(i),g,k,H,h(j),l);
a1(i,j)=acos(x(1))*180/pi;
end
end
a1即为出手角度。如取初值x0[0.7，1]，得到。计算结果见表4-1-4（将不考虑阻力的结果重写在最后两列，以作比较）。
表4-1-4 对不同出手速度和出手高度的出手角度
v(m/s）
h (m)
(度) (度)
(度) (度)
8.0
1.8
1.9
2.0
2.1
60.7869 43.5424
61.6100 41.5693
62.3017 39.7156
62.9012 37.9433
62.4099 42.7925
63.1147 40.9188
63.7281 39.1300
64.2670 37.4017
8.5
1.8
1.9
2.0
2.1
66.5719 37.7905
66.9244 36.2870
67.2505 34.7982
67.5541 33.3209
67.6975 37.5049
68.0288 36.0075
68.3367 34.5214
68.6244 33.0444
9.0
1.8
1.9
2.0
2.1
70.1198 34.2736
70.3328 32.9087
70.5352 31.5428
70.7279 30.1756
71.0697 34.1327
71.2749 3207614
71.4700 31.3881
71.6561 30.0127
思考 如何计算入射角，并由此判断出手角度是否应该舍去。
§1．7 结果分析
（1）最小出手速度和出手角度（表4-1-1）
对应于最小出手速度是最小出手角度，他们均随着出手高度的增加而略有减少；出手速度一般不要小于8米/秒。
（2）出手速度和出手高度对出手角度的影响（表4-1-2）
速度一定时，出手高度越大时，出手角度越大，但是随着速度的增加，高度对角度的影响变小，这种影响在1左右；出手高度一定时，速度越大，出手角度也应越大，速度的影响在。
（3）出手角度和出手速度的允许偏差（表4-1-3）
总的看来，允许偏差都相当小。进一步分析可知，出手高度一定，速度越大，角度的允许偏差越小，而速度的允许偏差越大，且对角度的要求比对速度的要求严格；出手速度一定，高度越大，虽然也是角度的允许偏差越小，速度的允许偏差越大，但这时对角度和速度的要求都相对较低。
（4）空气阻力的影响对同样的出手速度和高度，不考虑阻力影响时出手角度要大一些（）。
还可将空气阻力的影响与前面讨论的命中偏前（后）的允许偏差作一对比。由（4-1-17）式可知阻力对的影响因子，因为=0.05（1/s），可由fsolve的x（2）得到），所以阻力对命中偏前（后）的影响不超过3%；而由（4-1-9）式允许偏差约为0.09（m）（以计），与=4.6（m ）的相对偏差约2%，二者是相当的。
思考 若出手高度和角度固定，考察阻力对出手速度的影响。
§2水塔流量估计
§2.1问题的提出某居民区有一供居民用水的圆柱形水塔，一般可以通过测量其水位来估计水的流量.但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量.通常水泵每天供水一两次，每次约两小时。
水塔是一个高12.2、直径17.4米的正圆柱.按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作。
表4-2-1是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的流水量，及一天的总用水量。
表4-2-1 水位测量记录（符号//表示水泵启动）
时刻（h）
水位(cm)
0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97
968 948 931 913 898 881 869 852 839 822
时刻（h）
水位(cm)
9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93
// // 1082 1050 1021 994 965 941 918 892
时刻（h）
水位(cm)
19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91
866 843 822 // // 1059 1035 1018
§2.2问题的分析流量是单位时间流出的水的体积。由于水塔为正圆柱形，横截面积是常数，在水泵不工作的时段，流量很容易从水位对时间的变化率算出，问题是如何估计水泵供水时段的流量。
水泵供水时段的流量只能靠供水时段前后的流量拟合得到。作为用于拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。这些流量大体上可由两种方法计算，一是直接对表4-2-1中的水位用数值微分算出各时段的流量，用他们拟合其他时刻或连续时间的流量；二是先用表中的数据拟合水位~时间函数，求导数即可得到连续时间的流量。一般来说数值微分的精度不高，何况测量记录还是不等距的，数值微分的计算尤为麻烦。下面我们用第二种方法处理。
有了任何时刻的流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，如由表4-2-1可知从到（小时）水位下降了968-822=146（厘米），乘以水塔的截面积就是这一时段的用水量。这个数值可以用来检验拟合的结果。
§2.3 模型假设
（1）流量只取决于水位差，与水位本身无关。按照Torricelli定律从小孔流出的流体的流速正比于水面高度的平方根，题目给出水塔的最高和最低水位分别为10.8米和8.2
米（设出口的水位为0），因为，可忽略水位对流速的影响。
（2）水泵第1次供水时段为到（小时），第2次供水时段为到（小时）。这是根据最高和最低水位分别为10.8米和8.2米，及表1的水位测量记录做出的假设，其中前3个时刻直接取自实测数据（精确到0.1小时），最后1个时刻来自每次供水约两小时的已知条件（从记录看，第2供水时段应在有记录的22.96小时之后不久结束）。
（3）水泵工作时单位时间的供水量大致为常数，这个常数大于单位时间的平均流量
（4）流量是（对时间的）连续函数。
（5）流量与水泵是否工作无关。
（6）由于水塔截面积是常数，为简单起见，计算中将流量定义作单位时间流出的水的高度，即水位对时间变化率的绝对值（水位是下降的），最后给出结果时再乘以即可。
§2.4 流量估计
（1）拟合水位~时间函数从表4-2-1测量记录看，一天有两个水泵供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段（以下称第1时段到，第2时段到，和第3时段以后）。对第1、2时段的测量数据直接分别做多项式函数拟合，得到水位函数。为使拟合曲线比较光滑，多项式次数不要太高，一般在3~6。由于第3时段只有3个测量记录，无法对这一时段的水位做出较好的拟合。
（2）确定流量~时间函数对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内。
（3）一天总用水量的估计总用水量等于两个水泵不工作时段和两个供水时段用水量之和，他们都可以由流量对时间的积分得到。
§2.5 算法设计与编程
（1）拟合第1、2时段的水位，并导出流量设为已输入的时刻和水位测量记录（水泵启动的4个时刻不输入），第1时段各时刻的流量可如下得到：
c1=ployfit（t（1：10），h（1：10），3）；%用3次多项式拟合第1时段水位，c1输出3 次多项式的系数
a1=ployder（c1）；%a1输出多项式（系数为c1）导数的系数
tp1=0：0.1：9；
x1=-ployval（a1，tp1）；%x1输出多项式（系数为a1）在tp1点的函数值（取负后变为正值），即tp1时刻的流量类似地计算第2时段各时刻的流量。
（2）拟合水泵工作时段的流量在第1供水时段（）之前（即第1时段）和之后（即第2时段）各取几点，其流量已经得到，用他们拟合第1供水时段的流量。为使流量函数在和连续，我们简单地只取4个点，拟合3次多项式（即曲线必过这4个点），实现如下：
xx1=-polyval（a1，[8 9]）；%取第1时段在t=8，9的流量
xx2=-polyval（a2，[11 12]）；%取第2时段在t=11，12的流量
xx12=[xx1 xx2]；
c12=polyfit（[8 9 1 1 1 2]，xx12，3）；%拟合3次多项式
tp12=9：0.1：11；
x12=polyval（c12，tp12）；%x12输出第1供水时段各时刻的流量在第2供水时段之前取两点的流量，在该时段之后（第3时段）仅有3个水位记录，我们用差分得到流量，然后用这4个数值拟合第2供水时段的流量如下：
dt3=diff（t（26：28））；%用最后3个时刻的两两之差
dh3=diff（h（26：28））；%用最后3个时刻的两两之差
dht3=[-dh3./dt3；%t（26）和t（27）的流量
t3=[20 20.8 t（26） t（27）]；
xx3=[-polyval（a2，t3（1：2）），dht3]；%取t3各时刻的流量
c3=polyfit（t3，xx3，3）；%拟合3次多项式
tp3=20.8：0.1：24；
x3=polyval（c3，tp3）；%x3输出第2供水时段（外推至t=24）各时刻的流量
（3）一天总用水量的估计第1、2时段和第1、2供水时段流量的积分之和，就是一天总用水量。虽然诸时段的流量已表为多项式函数，积分可以解析的算出，这里仍用数值积分计算如下：
y1=0.1*trapz（x1）；%第1时段用水量（仍按高度计），0.1为积分步长
y2=0.1*trapz（x2）；%第2时段用水量
y12=0.1*trapz（x12）；%第1供水时段用水量
y3=0.1*trapz（x3）；%第2供水时段用水量
y=（y1+ y2+ y12+ y3）*237.8*0.01%一天总用水量（）
（4）流量及总用水量的检验计算出的各时刻的流量可用水位记录的数值微分来检验。用水量y1可用第1时段的水位测量记录中的下降高度968-822=146来检验，类似地，y2用1082-822=260检验。
供水时段流量的一种检验办法如下：供水时段的用水量加上水位升高值260是该时段水泵注入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供水时段水泵的功率应大致相等。第1、2时段水泵的功率可计算如下：
p1=（y12+260）/2；%第1供水时段水泵的功率（水量仍以高度计）
tp2=20.8：1：23；
xp2=ployval（c3，tp2）；%xp2输出第2供水时段各时刻的流量
p2=（0.1*trapz（xp2）+260）/2.2；第2供水时段水泵的功率（水量仍以高度计）
§2.6 计算结果从上面的分析和算法设计过程看，计算结果与各时段所用拟合的多项式的次数有关，下面只给出拟合第1、2时段水位函数时，采用不同次数的多项式所得流量及总用水量的结果。表4-1-2是各时段的用水量，一天总用水量及两个供水时段水泵的功率。
表4-1-22 各时段和一天总用水量 及两个供水时段水泵的功率（符号和单位见上）
(n1,n2)
y1 y2 y12 y3
y
p1 p2
(3,4)
146.18 258.10 48.50 78.50
1263.4
154.25 143.59
(5,6)
146.52 257.76 46.13 76.30
1252.5
153.07 142.67
§2.7 分析及改进由表4-1-2可以看出，第1时段用水量与水位测量记录中的下降高度146相差无几，第2时段用水量与记录中的下降高度260相差无几，所以数据拟合、数值积分的精度是足够的。对不同次数的拟合多项式第1、2供水时段用水量相差稍大，两供水时段水泵的功率亦有差别，都说明供水时段用3次曲线通过4点的做法不够好，应该多取几点拟合，但要注意让流量曲线在不同时段相接处保持连续。
由图可以看出，流量曲线与原始记录基本上相吻合，零点到十点钟流量很低，十点到下午三点是用水高峰，全天流量平均在22（cm）/h）左右。若按这个平均流量计算，一天总用水量应为，与表2中的结果非常接近。
§3 钢管订购和运输
§3.1问题的提出
要铺设一条的输送天然气的主管道,如图一所示(见下页)。经筛选后可以生产这种主管道钢管的钢厂有。图中粗线表示铁路，单细线表示公路，双细线表示要铺设的管道(假设沿管道或者原来有公路，或者建有施工公路)，圆圈表示火车站，每段铁路、公路和管道旁的阿拉伯数字表示里程(单位km)。
为方便计，1km主管道钢管称为1单位钢管。
一个钢厂如果承担制造这种钢管，至少需要生产500个单位。钢厂在指定期限内能生产该钢管的最大数量为个单位，钢管出厂销价1单位钢管为万元，如下表：

1
2
3
4
5
6
7

800
800
1000
2000
2000
2000
3000

160
155
155
160
155
150
160
1单位钢管的铁路运价如下表：
里程(km)
≤300
301～350
351～400
401～450
451～500
运价(万元)
20
23
26
29
32
里程(km)
501～600
601～700
701～800
801～900
901～1000
运价(万元)
37
44
50
55
60
1000km以上每增加1至100km运价增加5万元。
公路运输费用为1单位钢管每公里0.1万元（不足整公里部分按整公里计算）。
钢管可由铁路、公路运往铺设地点（不只是运到点，而是管道全线）。
（1）请制定一个主管道钢管的订购和运输计划，使总费用最小（给出总费用)。
（2）请就（1）的模型分析：哪个钢厂钢管的销价的变化对购运计划和总费用影响最大，哪个钢厂钢管的产量的上限的变化对购运计划和总费用的影响最大，并给出相应的数字结果。
（3）如果要铺设的管道不是一条线，而是一个树形图，铁路、公路和管道构成网络，请就这种更一般的情形给出一种解决办法，并对图二按（1）的要求给出模型和结果。
§3.2模型的假设:
（1）铺设钢管各线段运输费以1公里为单位,在1公里内运输费不计.
（2）钢管的铺设情况以从某一头向另一头铺设为准,某一段两头对铺设的情况不考虑.
（3）在钢管厂承担生产任务的情况下,其供应量至少为500个单位
§3.2符号约定
 第段长度；
 第段内的运费；
 第段内用第个工厂单位钢管的最小费用;
 第段内用第个工厂钢管的数量；
 第个工厂单位钢管的生产价格；
第个工厂生产钢管的上限；
§3.3问题的分析在问题(1)中,铺设钢管的总费用包括三个部分,(一)钢管的生产成本 (二) 钢管运到各段的费用 (三) 钢管铺设时的运输费用,根据模型假设2) 可以得知钢管铺设时的运输费用是一定值,并可用等差数列求出 ; 而钢管的生产成本与钢管运到各段的费用可以综合起来考虑,将单位生产价格与单位运费（以下简称运费） 相加得出单位钢管运到各段的费用。由此，问题归结为求各厂单位钢管运到各段的最小费用及求各厂运到各段的最优运量问题。
在问题（2）中可将各厂的最大生产上限和生产价格逐一进行变动，在一厂变动时，其他厂不变，这样可以看出各厂的生产上限与生产价格变动对总费用的影响。
在问题（3）中，处理问题的方法与问题（1）类似.
§3.4模型的建立与处理我们的目标函数是总费用，它包含两项：从工厂到各段的运费及其材料费之和和铺设时的运费,那么,根据假设易知应为一定值，即(其中，0.1为单位钢管的每公里运费)。对于,设单位生产价格与单位运费相加之和为单位运费,则；所以
 (4-3-1)
此目标函数的限定条件为：
1、各厂运往各段的钢管量之和为此段铺设钢管的总量:即
 （4-3-2）
2、各路段所需各厂的钢管量之和应小于此厂的最大生产量，大于其最小生产量，即
 （4-3-3）
§3.4.1问题一
（1）的计算：
根据公式，易求得：
 （4-3-4）
（2）的计算：
对于有多种计算方法：观察法、枚举法和动态规划法。
由于此题中给出的运输路线较少易于观察，比动态规划方法简单，所以在此我们采用了观察法和枚举法。通过对路段不同路径的运输费用进行计算比较，最终得出工厂到各路段的单位产品的最优费用，得出的矩阵为一的矩阵，其中行下标表示第路段，列下标表示第厂，具体矩阵如下：

从而可以利用线性规划中的单纯形法可求知，使达到最小，由于运算复杂，在此就不做考虑用MATLAB软件中的有关求解线性规划的方法求得，首先考虑各厂均供应，即，利用程序hyw.01可解得：
 （4-3-5）
具体各厂供应量如下：
单位







供应量
800
800
1000
500
854
717
500
在运算过程中可得到生产量全部供应给了这一段，介于实际情况考虑（即运往的总费用更小，且供应量充足），故可省去不考虑，这样使全局更优，即，用程序hyw.02计算可得：

具体各厂供应量如下：
单位







供应量
800
800
1000
500
767
1304
0
同样可知供应给了这两段，然而的材料价格及运输情况相对更高（且从上面数据可以发现即使使用代替其供应量不会有影响），从而在的情况下再次考虑，这样会使全局最优，用程序hyw.03计算可得：

具体各厂供应量如下：
单位







供应量
800
800
1000
0
1201
1370
0
各厂供应量（即（）代表路段，代表供应厂）分配如下：
供应厂路段







A1----A2
0
104
0
0
0
0
0
A2----A3
0
125.7778
67.4444
0
107.7778
0
0
A3----A4
272.5
184.6111
126.2778
0
166.6111
0
0
A4----A5
128.5
184.6111
126.2778
0
166.6111
0
0
A5----A6
194
0
0
0
0
0
0
A6----A7
205
0
0
0
0
0
0
A7----A8
0
201
0
0
0
0
0
A8----A9
0
0
680
0
0
0
0
A9----A10
0
0
0
0
240
240
0
A10---A11
0
0
0
0
300
0
0
A11---A12
0
0
0
0
220
0
0
A12---A13
0
0
0
0
0
210
0
A13---A14
0
0
0
0
0
420
0
A14---A15
0
0
0
0
0
500
0
以上计算及统计情况见程序hyw.03：故由（4）（7）可得：
min(W)=min()==1170626.1+122387.2=1293013.3
§3.4.2问题二分别将各厂的生产价格与生产上限变化一个单位，看他们对总费用的影响。具体可以用销价和生产上限变化一个单位，从而得到影响大小。
见下表：
S变化为1时








103
35
30
0
0
0
0
P变化为1时








800
800
1000
500
854
717
500
上表为销价变化及生产上限变化对总费用的影响情况，可以看出厂的产量上限变化对购运计划总费用影响最大,厂的销价变化对总费用的影响最大。
§3.4.3 问题三对于更一般的情况，如树形图，考虑情况可与（1）类似。 在厂的取舍过程中亦同。具体模型与求解的最佳方案一致，从而再次重复的求解过程，用程序hyw.04可得：=1296491.1;易知=128973.8;从而:
min(W)=min()+=1296491.1+128973.8=1425464.9
具体各厂供应总量如下：
单位







供应量
800
800
1000
0
1673
1630
0
各厂供应量（即（）i代表路段，j代表供应厂）分配如下：
供应厂路段







A1----A2
0
104
0
0
0
0
0
A2----A3
0
125.7778
53.4444
0
121.7778
0
0
A3----A4
272.5
184.6111
22.2778
0
180.6111
0
0
A4----A5
128.5
184.6111
22.2778
0
180.6111
0
0
A5----A6
194
0
0
0
0
0
0
A6----A7
205
0
0
0
0
0
0
A7----A8
0
201
0
0
0
0
0
A8----A9
0
0
680
0
0
0
0
A9----A10
0
0
0
0
480
0
0
A10---A11
0
0
0
0
300
0
0
A11---A12
0
0
0
0
220
0
0
A12---A13
0
0
0
0
0
210
0
A13---A14
0
0
0
0
0
420
0
A14---A15
0
0
0
0
0
500
0
A15---A16
0
0
42
0
0
0
0
A16---A17
0
0
0
0
0
10
0
A17---A18
0
0
0
0
0
130
0
A18---A19
0
0
0
0
190
0
0
A19---A20
0
0
0
0
0
260
0
A20---A21
0
0
0
0
0
100
0
§3.5 模型的优点与改进优点：本模型适用范围较广，在辅助软件的求解下，结果易得出，管道数量的变化对算法几乎没有影响。
改进：在实际情况中，可考虑从某一段从两头对铺，使模型更优化,对于供应量较多，运费和单价较低情况下应就近考虑(如：) 。在无软件帮助的情况下，可考虑一些不可能的情况，使计算规模缩小。
附录：
hw01:
clear;
s=[800 800 1000 2000 2000 2000 3000];
p1=[160 155 150 160 155 150 160];
q=[160.3 140.2 98.6 38 20.5 3.1 3.1 21.2 64.2 92 96 106 121.1 128;
205.3 190.2 171.6 111 95.5 86 71.2 71.2 114.2 142 146 156 171.1 178;
220.3 200.2 181.6 121 105.5 96 86.2 48.2 48.2 82 86 96 111.2 118;
250.3 235.2 216.6 156 130.5 133.1 116.2 84.2 62 51 51 61 76.2 83;
245.3 225.2 206.6 146 140.5 121 111.2 79.2 57 33 33 51 71.2 73;
255.3 235.2 216.6 156 140.5 131 121.1 84.2 62 51 45 26.2 11 11;
275.3 245.2 226.6 166 150.5 141 131.2 99.2 77 66 56 38.2 26 2];
p=p1'*ones(1,14);
q=q+p;
c=[];
for i=1:14
c=[c q(:,i)'];
end
d=zeros(14,98);
for i=1:14
for j=7*i-6:7*i
d(i,j)=1;
end
end
B=[];
for i=1:14
B=[B eye(7)];
i=i+1;
end
A=[d;B;-B];
b1=[104,301,750,606,194,205,201,680,480,300,220,210,420,500];
b2=[s(1),s(2),s(3),s(4),s(5),s(6),s(7)];
b3=-500*ones(1,7);
b=[b1';b2';b3'];
v1=zeros(98,1);
x=lp(c,A,b,v1,[],[],14)
M=vpa(c*x,8)
S=B*x
hw02:
clear;
s=[800 800 1000 2000 2000 2000];
p1=[160 155 150 160 155 150];
q=[160.3 140.2 98.6 38 20.5 3.1 3.1 21.2 64.2 92 96 106 121.1 128;
205.3 190.2 171.6 111 95.5 86 71.2 71.2 114.2 142 146 156 171.1 178;
220.3 200.2 181.6 121 105.5 96 86.2 48.2 48.2 82 86 96 111.2 118;
250.3 235.2 216.6 156 130.5 133.1 116.2 84.2 62 51 51 61 76.2 83;
245.3 225.2 206.6 146 140.5 121 111.2 79.2 57 33 33 51 71.2 73;
255.3 235.2 216.6 156 140.5 131 121.1 84.2 62 51 45 26.2 11 11];
p=p1'*ones(1,14);
q=q+p;
c=[];
for i=1:14
c=[c q(:,i)'];
end
d=zeros(14,84);
for i=1:14
for j=6*i-5:6*i
d(i,j)=1;
end
end
B=[];
for i=1:14
B=[B eye(6)];
i=i+1;
end
A=[d;B;-B];
b1=[104,301,750,606,194,205,201,680,480,300,220,210,420,500];
b2=[s(1),s(2),s(3),s(4),s(5),s(6)];
b3=-500*ones(1,6);
b=[b1';b2';b3'];
v1=zeros(84,1);
x=lp(c,A,b,v1,[],[],14)
M=vpa(c*x,8)
S=B*x
hw03:
clear;
s=[800 800 1000 2000 2000];
p1=[160 155 150 155 150];
q=[160.3 140.2 98.6 38 20.5 3.1 3.1 21.2 64.2 92 96 106 121.1 128;
205.3 190.2 171.6 111 95.5 86 71.2 71.2 114.2 142 146 156 171.1 178;
220.3 200.2 181.6 121 105.5 96 86.2 48.2 48.2 82 86 96 111.2 118;
245.3 225.2 206.6 146 140.5 121 111.2 79.2 57 33 33 51 71.2 73;
255.3 235.2 216.6 156 140.5 131 121.1 84.2 62 51 45 26.2 11 11];
p=p1'*ones(1,14);
q=q+p;
c=[];
for i=1:14
c=[c q(:,i)'];
end
d=zeros(14,70);
for i=1:14
for j=5*i-4:5*i
d(i,j)=1;
end
end
B=[];
for i=1:14
B=[B eye(5)];
i=i+1;
end
A=[d;B;-B];
b1=[104,301,750,606,194,205,201,680,480,300,220,210,420,500];
b2=[s(1),s(2),s(3),s(4),s(5)];
b3=-500*ones(1,5);
b=[b1';b2';b3'];
v1=zeros(70,1);
x=lp(c,A,b,v1,[],[],14)
M=vpa(c*x,8)
S=B*x
hw04:
clear;
s=[800 800 1000 2000 2000];
p1=[160 155 150 155 150];
q1=[160.3 140.2 98.6 38 20.5 3.1 3.1 21.2 64.2 92 96 106 121.1 128;
205.3 190.2 171.6 111 95.5 86 71.2 71.2 114.2 142 146 156 171.1 178;
220.3 200.2 181.6 121 105.5 96 86.2 48.2 48.2 82 86 96 111.2 118;
245.3 225.2 206.6 146 140.5 121 111.2 79.2 57 33 33 51 71.2 73;
255.3 235.2 216.6 156 140.5 131 121.1 84.2 62 51 45 26.2 11 11];
q2=[60 95 95 95 105 115;
110 145 145 145 155 165;
44 85 85 85 65 105;
75 32 32 32 50 65;
80 32 36 50 10 0];
q=[q1 q2];
p=p1'*ones(1,20);
q=q+p;
c=[];
for i=1:20 c=[c q(:,i)'];
end
d=zeros(20,100);
for i=1:20
for j=5*i-4:5*i
d(i,j)=1;
end
end
B=[];
for i=1:20
B=[B eye(5)];
i=i+1;
end
A=[d;B;-B];
b1=[104,301,750,606,194,205,201,680,480,300,220,210,420,500,42,10,130,190,260,100];
b2=[s(1),s(2),s(3),s(4),s(5)];
b3=-500*ones(1,5);
b=[b1';b2';b3'];
v1=zeros(100,1);
x=lp(c,A,b,v1,[],[],20)
M=vpa(c*x,8)
S=B*x
大作业：某服务部门一周中每天需要不同数目的雇员：周一至周四每天至少需要50人，周五和周日每天至少80人，周六至少90人，现规定应聘者需连续工作5天
1）、试建立模型确定聘用方案：即周一到周日每天聘用多少人，使在满足需要的条件下聘用总人数最少？
2）、上面指的是全时雇员（8小时/天）如果可以用两个临时聘用的半时雇员（4小时/天）代替一个全时雇员，但规定半时雇员的工作量不得超过总工作量的，全时雇员和半时雇员每小时酬金分别为5元和3元，试建立模型确定聘用方案：在满足需要的条件下所付酬金总额最小。