简单的数学建模

数学建模1

数学模型（MathematicalModel）是用数学符号、数学式子、程序、图形等对实际课题本质属性的抽象而又简洁的刻划，它或能解释某些客观现象，或能预测未来的发展规律，或能为控制某一现象的发展提供某种意义下的最优策略或较好策略。

数学建模（MathematicalModeling）

应用知识从实际课题中抽象、提炼出数学模型的过程。§1.1

数学模型与数学建模2

1.了解问题的实际背景，明确建模目的，收集掌握必要的数据资料。

2.在明确建模目的，掌握必要资料的基础上，通过对资料的分析计算，找出起主要作用的因素，经必要的精炼、简化，提出若干符合客观实际的假设。

3.在所作假设的基础上，利用适当的数学工具去刻划各变量之间的关系，建立相应的数学结构——即建立数学模型。

4.模型求解。

5.模型的分析与检验。

在难以得出解析解时，也应当借助计算机求出数值解。

§1.2

数学建模的一般步骤实体信息(数据)假设建模求解验证应用3§1.3

数学模型的分类分类标准具体类别对某个实际问题了解的深入程度白箱模型、灰箱模型、黑箱模型模型中变量的特征连续型模型、离散型模型或确定性模型、随机型模型等建模中所用的数学方法初等模型、微分方程模型、差分方程模型、优化模型等研究课题的实际范畴人口模型、生态系统模型、交通流模型、经济模型、基因模型等4例1

某人平时下班总是按预定时间到达某处，然然后他妻子开车接他回家。有一天，他比平时提早了三十分钟到达该处，于是此人就沿着妻子来接他的方向步行回去并在途中遇到了妻子，这一天，他比平时提前了十分钟到家，问此人共步行了多长时间？

§1.4一些简单实例

似乎条件不够哦。。

换一种想法，问题就迎刃而解了。假如他的妻子遇到他后仍载着他开往会合地点，那么这一天他就不会提前回家了。提前的十分钟时间从何而来？

显然是由于节省了从相遇点到会合点，又从会合点返回相遇点这一段路的缘故，故由相遇点到会合点需开5分钟。而此人提前了三十分钟到达会合点，故相遇时他已步行了二十五分钟。

请思考一下，本题解答中隐含了哪些假设？5例2

交通灯在绿灯转换成红灯时，有一个过渡状态——亮一段时间的黄灯。请分析黄灯应当亮多久。设想一下黄灯的作用是什么，不难看出，黄灯起的是警告的作用，意思是马上要转红灯了，假如你能停住，请立即停车。停车是需要时间的，在这段时间内，车辆仍将向前行驶一段距离L。这就是说，在离街口距离为L处存在着一条停车线（尽管它没被画在地上），见图。对于那些黄灯亮时已过线的车辆，则应当保证它们仍能穿过马路。

马路的宽度D是容易测得的，问题的关键在于L的确定。为确定L，还应当将L划分为两段：L1和L2，其中L1是司机在发现黄灯亮及判断应当刹车的反应时间内驶过的路程，L2为刹车制动后车辆驶过的路程。L1较容易计算，交通部门对司机的平均反应时间t1早有测算，反应时间过长将考不出驾照），而此街道的行驶速度v也是交管部门早已定好的，目的是使交通流量最大，可另建模型研究，从而L1=v*t1。刹车距离L2既可用曲线拟合方法得出，也可利用牛顿第二定律计算出来。黄灯究竟应当亮多久现在已经变得清楚多了。第一步，先计算出L应多大才能使看见黄灯的司机停得住车。第二步，黄灯亮的时间应当让已过线的车顺利穿过马路，即T至少应当达到（L+D）/v。

DL6

2.初等模型举例7

常见类型定性模型经验公式（拟合、插值）量纲分析比例模型8§2.1

崖高的估算假如你站在崖顶且身上带着一只具有跑表功能的计算器，你也许会出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度，假定你能准确地测定时间，你又怎样来推算山崖的高度呢，请你分析一下这一问题。我有一只具有跑表功能的计算器。9方法一假定空气阻力不计，可以直接利用自由落体运动的公式来计算。例如，设t=4秒，g=9.81米/秒2，则可求得h≈78.5米。

我学过微积分，我可以做得更好，呵呵。

10除去地球吸引力外，对石块下落影响最大的当属空气阻力。根据流体力学知识，此时可设空气阻力正比于石块下落的速度，阻力系数K为常数，因而，由牛顿第二定律可得：

令k=K/m，解得

代入初始条件v(0)=0，得c=－g/k，故有

再积分一次，得：

11若设k=0.05并仍设t=4秒，则可求得h≈73.6米。

听到回声再按跑表，计算得到的时间中包含了反应时间

进一步深入考虑不妨设平均反应时间为0.1秒，假如仍设t=4秒，扣除反应时间后应为3.9秒，代入式①，求得h≈69.9米。

①多测几次，取平均值再一步深入考虑代入初始条件h(0)=0，得到计算山崖高度的公式：

将e-kt用泰勒公式展开并令k→0+

，即可得出前面不考虑空气阻力时的结果。12还应考虑回声传回来所需要的时间。为此，令石块下落的真正时间为t1，声音传回来的时间记为t2，还得解一个方程组：这一方程组是非线性的，求解不太容易，为了估算崖高竟要去解一个非线性主程组似乎不合情理

相对于石块速度，声音速度要快得多，我们可用方法二先求一次

h，令t2=h/340，校正t，求石块下落时间t1≈t-t2将t1代入式①再算一次，得出崖高的近似值。例如，若h=69.9米，则t2≈0.21秒，故t1≈3.69秒，求得h≈62.3米。13§2.2

录像带还能录多长时间录像机上有一个四位计数器，一盘180分钟的录像带在开始计数时为0000，到结束时计数为1849，实际走时为185分20秒。我们从0084观察到0147共用时间3分21秒。若录像机目前的计数为1428，问是否还能录下一个60分钟的节目？14rθRl由得到又因和得

积分得到即从而有我们希望建立一个录像带已录像时间t与计数器计数n之间的函数关系。为建立一个正确的模型，首先必须搞清哪些量是常量，哪些量是变量。首先，录像带的厚度W是常量，它被绕在一个半径为r的园盘上，见图。磁带转动中线速度v显然也是常数，否则图象声音必然会失真。此外，计数器的读数n与转过的圈数有关，从而与转过的角度θ成正比。15rθRl

此式中的三个参数W、v和r均不易精确测得，虽然我们可以从上式解出t与n的函数关系，但效果不佳，故令则可将上式简化为：故令上式又可化简记成t=an2+bn

16t=an2+bn

rθRl上式以a、b为参数显然是一个十分明智的做法，它为公式的最终确立即参数求解提供了方便。将已知条件代入，得方程组：从后两式中消去t1，解得a=0.0000291,b=0.04646,故t=0.0000291n2+0.04646n，令n=1428，得到t=125.69（分）由于一盒录像带实际可录像时间为185.33分，故尚可录像时间为59.64分，已不能再录下一个60分钟的节目了。

17在解决实际问题时，注意观察和善于想象是十分重要的，观察与想象不仅能发现问题隐含的某些属性，有时还能顺理成章地找到解决实际问题的钥匙。本节的几个例子说明，猜测也是一种想象力。没有合理而又大胆的猜测，很难做出具有创新性的结果。开普勒的三大定律（尤其是后两条）并非一眼就能看出的，它们隐含在行星运动的轨迹之中，隐含在第谷记录下来的一大堆数据之中。历史上这样的例子实在太多了。在获得了一定数量的资料数据后，人们常常会先去猜测某些结果，然后试图去证明它。猜测一经证明就成了定理，而定理一旦插上想象的翅膀，又常常会被推广出许多更为广泛的结果。即使猜测被证明是错误的，结果也决不是一无所获的失败而常常是对问题的更为深入的了解。§2.3最短路径与最速方案问题18例5（最短路径问题）

设有一个半径为r的圆形湖，圆心为O。A、B

位于湖的两侧，AB连线过O，见图。现拟从A点步行到B点，在不得进入湖中的限制下，问怎样的路径最近。

ABOr将湖想象成凸出地面的木桩，在AB间拉一根软线，当线被拉紧时将得到最短路径。根据这样的想象，猜测可以如下得到最短路径：过A作圆的切线切圆于E，过B作圆的切线切圆于F。最短路径为由线段AE、弧EF和线段FB连接而成的连续曲线（根据对称性，AE′，弧E′F′，F′B连接而成的连续曲线也是）。EFE′F′19以上只是一种猜测，现在来证明这一猜测是正确的。为此，先介绍一下凸集与凸集的性质。定义2.1（凸集）称集合R为凸集，若x1、x2∈R及λ∈[0，1]，总有λx1+（1+λ）x2∈R。即若x1、x2∈R，则x1、x2的连线必整个地落在R中。定理2.2（分离定理）对平面中的凸集R与R外的一点K，存在直线l,l

分离R与K，即R与K分别位于l的两侧（注：对一般的凸集R与R外的一点K，则存在超平面分离R与K），见图。klR下面证明猜想20猜测证明如下：（方法一）显然，由AE、EF、FB及AE′，E′F′，F′B围成的区域R是一凸集。利用分离定理易证最短径不可能经过R外的点，若不然，设Γ为最短路径，Γ过R外的一点M，则必存在直线l分离M与R，由于路径Γ是连续曲线，由A沿Γ到M，必交l于M1，由M沿Γ到B又必交l于M2。这样，直线段M1M2的长度必小于路径M1MM2的长度，与Γ是A到B的最短路径矛盾，至此，我们已证明最短路径必在凸集R内。不妨设路径经湖的上方到达B点，则弧EF必在路径F上，又直线段AE是由A至E的最短路径，直线FB是由F到B的最短路径，猜测得证。ABOrEFE′F′M1M2MΓl21还可用微积分方法求弧长，根据计算证明满足限止条件的其他连续曲线必具有更大的长度；此外，本猜测也可用平面几何知识加以证明等。根据猜测不难看出，例5中的条件可以大大放松，可以不必设AB过圆心，甚至可不必设湖是圆形的。例如对下图，我们可断定由A至B的最短路径必为l1与l2之一，其证明也不难类似给出。ABl1l2D到此为止，我们的研讨还只局限于平面之中，其实上述猜测可十分自然地推广到一般空间中去。1973年，J.W.Craggs证明了以上结果：若可行区域的边界是光滑曲面。则最短路径必由下列弧组成，它们或者是空间中的自然最短曲线，或者是可行区域的边界弧。而且，组成最短路径的各段弧在连接点处必定相切。22例6

一辆汽车停于A处并垂直于AB方向，此汽车可转的最小圆半径为R，求不倒车而由A到B的最短路径。解（情况1）若|AB|>2R，最短路径由弧AC与切线BC组成（见图①

）。（情况2）若|AB|<2R，则最短路径必居于图②（a）、（b）两曲线之中。可以证明，（b）中的曲线ABC更短。AR2RBRC①②ABoC（a）CABo1o2（b）23例7

驾驶一辆停于A处与AB成θ1角度的汽车到B处去，已知B处要求的停车方向必须与AB成θ2角，试找出最短路径（除可转的最小圆半径为R外，不受其他限止）。解根据Craggs定理并稍加分析可知，最短路径应在l1与l2中，见图，比较l1与l2的长度，即可得到最短路径。Al1l2Bθ2θ124最速方案问题例8

将一辆急待修理的汽车由静止开始沿一直线方向推至相隔S米的修车处，设阻力不计，推车人能使车得到的推力f满足:-B≤f≤A,f>0为推力，f<0为拉力。问怎样推车可使车最快停于修车处。

设该车的运动速度为υ=υ(t)，根据题意，υ(0)=υ(T)=0，其中T为推车所花的全部时间。由于-B≤f≤A，且f=mυ′，可知-b≤υ′≤a（其中m为汽车质量，a=A/m,b=B/m）。据此不难将本例归纳为如下的数学模型：

minT

υ(0)=υ(T)=025此问题为一泛函极值问题，求解十分困难，为得出一个最速方案。我们作如下猜测：猜测最速方案为以最大推力将车推到某处，然后以最大拉力拉之，使之恰好停于修车处，其中转换点应计算求出证明设υ=υ(t)为在最速推车方案下汽车的速度，则有。设此方案不同于我们的猜测。现从O点出发，作射线y=at；从（t,0）出发，作直线y=-b(t-T)交y=at于A，由于，曲线υ=υ(t)必位于三角形区域DAT的内部，从而有ΔOAT的面积大于S。在O到T之间任取一点T′

，过T′作AT的平行线交OA于A′

。显然ΔOA′T′的面积S（T′）是T′的连续函数，当T′=0时S（0）=0，当T′=T时，S（T）>S，故由连续函数的性质存在某T′<T，S（T′）=S但这一结果与υ=υ(t)是最优方案下的车速的假设矛盾，因为用我们猜测的推车方法推车，只需T′时间即可将车推到修车处，而T′<T。oυtAT′TA′Sy=aty=-b(t-T)26

3.微分方程建模27

常用技巧工程师原则房室系统建模竞争项的统计筹算律集中参数法与分布参数法灵敏度分析稳定性分析28§3.1

为什么要用三级火箭来发射人造卫星构造数学模型，以说明为什么不能用一级火箭而必须用多级火箭来发射人造卫星？为什么一般都采用三级火箭系统？1、为什么不能用一级火箭发射人造卫星?

（1）卫星能在轨道上运动的最低速度假设：（i）卫星轨道为过地球中心的某一平面上的圆，卫星在此轨道上作匀速圆周运动。（ii）地球是固定于空间中的均匀球体，其它星球对卫星的引力忽略不计。分析：根据牛顿第三定律，地球对卫星的引力为:在地面有:得:k=gR2

R为地球半径，约为6400公里故引力:假设(ii)29dmm-dmvu-v假设(i)卫星所受到的引力也就是它作匀速圆周运动的向心力故又有:从而:设g=9.81米/秒2，得:

卫星离地面高度

(公里)卫星速度

(公里/秒)10020040060080010007.807.697.587.477.377.86（2）火箭推进力及速度的分析假设：火箭重力及空气阻力均不计分析：记火箭在时刻t的质量和速度分别为m(t)和υ(t)有：记火箭喷出的气体相对于火箭的速度为u（常数），由动量守恒定理：υ0和m0一定的情况下，火箭速度υ(t)由喷发速度u及质量比决定。

故：由此解得：(3.11)

30（2）火箭推进力及速度的分析现将火箭——卫星系统的质量分成三部分：（i）mP（有效负载，如卫星）（ii）mF（燃料质量）（iii）mS（结构质量——如外壳、燃料容器及推进器）。最终质量为mP+mS

，初始速度为0，所以末速度：根据目前的技术条件和燃料性能，u只能达到3公里/秒，即使发射空壳火箭，其末速度也不超过6.6公里/秒。目前根本不可能用一级火箭发射人造卫星火箭推进力在加速整个火箭时，其实际效益越来越低。如果将结构质量在燃料燃烧过程中不断减少，那么末速度能达到要求吗？312、理想火箭模型假设：

记结构质量mS在mS+mF中占的比例为λ，假设火箭能随时抛弃无用的结构，结构质量与燃料质量以λ与（1-λ）的比例同时减少。建模:

由

得到：解得：

理想火箭与一级火箭最大的区别在于，当火箭燃料耗尽时，结构质量也逐渐抛尽，它的最终质量为mP，所以最终速度为：

只要m0足够大，我们可以使卫星达到我们希望它具有的任意速度。考虑到空气阻力和重力等因素，估计（按比例的粗略估计）发射卫星要使υ=10.5公里/秒才行，则可推算出m0/mp约为51,即发射一吨重的卫星大约需要50吨重的理想火箭323、理想过程的实际逼近——多级火箭卫星系统

记火箭级数为n，当第i级火箭的燃料烧尽时，第i+1级火箭立即自动点火，并抛弃已经无用的第i级火箭。用mi表示第i级火箭的质量，mP表示有效负载。先作如下假设：

（i）设各级火箭具有相同的λ,即i级火箭中λmi为结构质量，（1-λ）mi为燃料质量。

（ii）设燃烧级初始质量与其负载质量之比保持不变，并记比值为k。考虑二级火箭：

由3.11式，当第一级火箭燃烧完时，其末速度为：当第二级火箭燃尽时，末速度为：该假设有点强加的味道，先权作讨论的方便吧33又由假设（ii），m2=kmP，m1=k(m2+mP)，代入上式，仍设u=3公里/秒，且为了计算方便，近似取λ=0.1，则可得：要使υ2=10.5公里/秒，则应使:即k≈11.2，而:类似地，可以推算出三级火箭：

在同样假设下:

要使υ3=10.5公里/秒，则(k+1)/(0.1k+1)≈3.21，k≈3.25，而（m1+m2+m3+mP）/mP≈77。三级火箭比二级火箭几乎节省了一半是否三级火箭就是最省呢？最简单的方法就是对四级、五级等火箭进行讨论。34考虑N级火箭：

记n级火箭的总质量（包含有效负载mP）为m0

，在相同的假设下可以计算出相应的m0/mP的值，见表3-2n（级数）12345…∞（理想）

火箭质量（吨）/149776560…50表3-2由于工艺的复杂性及每节火箭都需配备一个推进器，所以使用四级或四级以上火箭是不合算的，三级火箭提供了一个最好的方案。当然若燃料的价钱很便宜而推进器的价钱很贵切且制作工艺非常复杂的话，也可选择二级火箭。354、火箭结构的优化设计3中已经能说过假设(ii)有点强加的味道；现去掉该假设，在各级火箭具有相同λ的粗糙假设下，来讨论火箭结构的最优设计。W1=m1+…+mn+mP

W2=m2+…+mn+mP……Wn=mn+mPWn+1=mP记应用（3.11）可求得末速度：记则又问题化为，在υn一定的条件下，求使k1k2…kn最小

解条件极值问题：或等价地求解无约束极值问题：可以解出最优结构设计应满足：火箭结构优化设计讨论中我们得到与假设（ii）相符的结果，这说明前面的讨论都是有效的！36§3.2

药物在体内的分布

何为房室系统？

在用微分方程研究实际问题时，人们常常采用一种叫“房室系统”的观点来考察问题。根据研究对象的特征或研究的不同精度要求，我们把研究对象看成一个整体（单房室系统）或将其剖分成若干个相互存在着某种联系的部分（多房室系统）。

房室具有以下特征：它由考察对象均匀分布而成，房室中考察对象的数量或浓度（密度）的变化率与外部环境有关，这种关系被称为“交换”且交换满足着总量守衡。在本节中，我们将用房室系统的方法来研究药物在体内的分布。在下一节中，我们将用多房室系统的方法来研究另一问题。交换环境内部单房室系统均匀分布37

药物的分解与排泄（输出）速率通常被认为是与药物当前的浓度成正比的，即：药物分布的单房室模型

单房室模型是最简单的模型，它假设：体内药物在任一时刻都是均匀分布的，设t时刻体内药物的总量为x(t)；系统处于一种动态平衡中，即成立着关系式：

药物的输入规律与给药的方式有关。下面，我们来研究一下在几种常见的给药方式下体内药体的变化规律。机体环境药物总量图3-8

假设药物均匀分布38情况1快速静脉注射机体环境只输出不输入房室其解为：药物的浓度：

与放射性物质类似，医学上将血浆药物浓度衰减一半所需的时间称为药物的血浆半衰期：负增长率的Malthus模型

在快速静脉注射时，总量为D的药物在瞬间被注入体内。设机体的体积为V，则我们可以近似地将系统看成初始总量为D，浓度为D/V，只输出不输入的房室，即系统可看成近似地满足微分方程：（3.12）

39情况2恒速静脉点滴机体环境恒定速率输入房室药物似恒速点滴方式进入体内，即:则体内药物总量满足：（x(0)=0）

（3.13）

这是一个一阶常系数线性方程，其解为：或易见：称为稳态血药浓度

对于多次点滴，设点滴时间为T1，两次点滴之间的间隔时间设为T2，则在第一次点滴结束时病人体内的药物浓度可由上式得出。其后T2时间内为情况1。故：（第一次）

0≤t≤T1

T1≤t≤T1

+T2

类似可讨论以后各次点滴时的情况，区别只在初值上的不同。第二次点滴起，患者体内的初始药物浓度不为零。40情况3口服药或肌注y(t)x(t)K1yK1x环境机体外部药物

口服药或肌肉注射时，药物的吸收方式与点滴时不同，药物虽然瞬间进入了体内，但它一般都集中与身体的某一部位，靠其表面与肌体接触而逐步被吸收。设药物被吸收的速率与存量药物的数量成正比，记比例系数为K1，即若记t时刻残留药物量为y(t)，则y满足：D为口服或肌注药物总量

因而：所以：解得：从而药物浓度：41

图3-9给出了上述三种情况下体内血药浓度的变化曲线。容易看出，快速静脉注射能使血药浓度立即达到峰值，常用于急救等紧急情况；口服、肌注与点滴也有一定的差异，主要表现在血药浓度的峰值出现在不同的时刻，血药的有效浓度保持时间也不尽相同。图3-9

我们已求得三种常见给药方式下的血药浓度C(t)，当然也容易求得血药浓度的峰值及出现峰值的时间，因而，也不难根据不同疾病的治疗要求找出最佳治疗方案。42

新药品、新疫苗在临床应用前必须经过较长时间的基础研究、小量试制、中间试验、专业机构评审及临床研究。当一种新药品、新疫苗研制出来后，研究人员必须用大量实验搞清它是否真的有用，如何使用才能发挥最大效用，提供给医生治病时参考。在实验中研究人员要测定模型中的各种参数，搞清血药浓度的变化规律，根据疾病的特点找出最佳治疗方案（包括给药方式、最佳剂量、给药间隔时间及给药次数等），这些研究与试验据估计最少也需要数年时间。在2003年春夏之交的SARS（非典）流行期内，有些人希望医药部门能赶快拿出一种能治疗SARS的良药或预防SARS的有效疫苗来，但这只能是一种空想。SARS的突如其来，形成了“外行不懂、内行陌生”的情况。国内权威机构一度曾认为这是“衣原体”引起的肺炎，可以用抗生素控制和治疗。但事实上，抗生素类药物对SARS的控制与治疗丝毫不起作用。以钟南山院士为首的广东省专家并不迷信权威，坚持认为SARS是病毒感染引起的肺炎，两个月后（4月16日），世界卫生组织正式确认SARS是冠状病毒的一个变种引起的非典型性肺炎（注：这种确认并非是由权威机构定义的，而是经对猩猩的多次实验证实的）。发现病原体尚且如此不易，要攻克难关，找到治疗、预防的办法当然就更困难了，企图几个月解决问题注定只能是一种不切实际的幻想。43

上述研究是将机体看成一个均匀分布的同质单元，故被称单房室模型，但机体事实上并不是这样。药物进入血液，通过血液循环药物被带到身体的各个部位，又通过交换进入各个器官。因此，要建立更接近实际情况的数学模型就必须正视机体部位之间的差异及相互之间的关联关系，这就需要多房室系统模型。IIIk12k21两房室系统图3-10

图3-10表示的是一种常见的两房室模型，其间的k12表示由室I渗透到室II的变化率前的系数，而k21则表示由室II返回室I的变化率前的系数，它们刻划了两室间的内在联系，其值应当用实验测定，使之尽可能地接近实际情况。当差异较大的部分较多时，可以类似建立多房室系统，即N房室系统44

4.离散优化模型45常用技巧计算复杂性分析算法设计

精确算法近似算法算法计算量估计、算法优劣比较46比较算法的好坏，从不同的角度出发，有各种不同的标准。在这里，我们仅就算法的计算速度作一个十分粗略的比较。

例1

（整理问题）给定n个实数a1,a2,…,an，要求将它整理成由小到大排列（或由大到小排列）的顺序：b1,b2,…,bn，b1≤b2≤…≤bn。（算法1）取出a1,a2,…,an中的最小者，令其为b1。从a1,a2,…,an中去除b1，在余下的n—1个数中选出最小者，令其为b2，…，直至得到b1,b2,…,bn。容易看出，为了排出b1,b2,…,bn，算法工作了次比较。（算法2）步0b1←a1

步1设已有b1,…,bk(1≤k<n)，将按两分法比较的方式把ak+1排入其中：若b1≤…≤bi≤ak+1≤bi+1≤…≤bk，令（b1,b2,…,bk,bk+1）←（b1,…,bi,ak+1,bi+1,…,bk）。若k+1<n，令k←k+1，返回步1。47我们来分析一下算法2的计算量：排出b1不必作比较，排出b2只需作一次比较，…，一般，在排ak+1时，设2r－1≤k＜2r，则只需作r次比较即可将ak+1安排在它应排的位置上。例如在排a8时，k=7，先和b4比，若a8>b4，可再和b6比（若a8<b4则和b2比），易见，只要比3次即可排入a8，由于r≤log2k+1，算法的总经较次数不超过令，，设计算机每秒可作C次比较，则算法1与算法2整理a1,a2,…,an所用的时间分别为

若n=100万，C=100万次/秒，则t1≈5.8天，而t2≈20秒。可见在较大规模的整理问题时，算法2明显优于算法1。

48现设一台每小时能作M次运算的计算机，并设求解某一问题已有了两个不同的算法。算法A对规模为n的实例约需作n2次运算而算法B则约需作2n次运算。于是，运用算法A在一小时内可解一个规模为的实例，而算法B则只能解一个规模为log2M的实例。两者的差别究竟有多大呢？让我们来对比一下。假如计算机每秒可作100万次运算，则算法A一小时大约可解一个规模为n=60000的实例，而算法B在一小时内只能解一个规模的实例且n每增加1，算法B求解时所化的时间就将增加1倍。例如算法B求解一个n=50的实例将连续运算357年多。现设计算速度提高了100倍，这对计算机来讲已是一个相当大的改进。此时算法A可解问题的规模增大了10倍，而算法B可解问题的规模只增加了log2100<7。前者可解问题的规模成倍成倍地增加而后者则几乎没有什么改变，今天无法求解的问题将来也很少有希望解决。

由于这一原因，人们对算法作了如下的分类：49（多项式算法）设A是求解某一问题D的一个算法，对D的一个规模为n的实例，用fA(D,n)表示用算法A求解这一实例所作的初等运算的次数。若存在一个多项式P(n)和一个正整数N，当n≥N时总有fA(D,n)≤P(n)（不论求解D的怎样的实例，只要其规模为n），则称算法A为求解问题D的一个多项式时间算法，简称多项式算法。定义4.1定义4.2（指数算法）设B是求解某一问题D的一个算法，fB(D,n)为用算法B求解D的一个规模为n的实例时所用的初等运算次数，若存在一个常数k>0，对任意正整数N，总可找到一个大于N的正整数n及D的一个规模为n的实例，对这一实例有fB(D,n)=O(2kn)，则称B为求解问题D的一个指数算法。多项式算法被称为是“好”的算法即所谓有效算法，而指数算法则一般被认为是“坏”算法，因为它只能求解规模极小的实例。50下面的表1列出了在规模大约为n时各类算法的计算量，而表8-15则反映了当计算机速度提高10倍、100倍时，各类算法在1小时内可求解的问题的模型的增长情况，（前三个是多项式时间的，后两个是指数时间的）表1几个多项式和指数时间复杂性函数增长情况算法要求的计算量规模n的近似值101001000n101001000nlogn336649966n31031061092n10241.27×10301.05×10301n!101584×10256751表21小时内可解的问题实例的规模算法要求的计算量用现在计算机用快10倍计算机用快100倍计算机nN110N1100N1nlognN28.2N267N2n3N32.15N34.64N32nN4N4+3.32N4+6.64n!N5≤N5+2≤N5+3由定义1知4n2与2n2都是O(n2)，nlogn+3n是O(nlogn)，我们在以后分析时间复杂性函数时也往往忽略常系数和增长速度较慢的项，因为前者可通过提高计算机速度来提高效率，而后者增长速度最快的项才是决定效率的关键因素。

六个最初的NP难问题52（1）（3—满足问题，简记3—SAT问题）每一个句子都是一个三项式的SAT问题，称为3—SAT问题。例如，就是3—SAT的一个实例。2）（三维匹配问题——3DM）X、Y、Z是三个不相交的集合，|X|=|Y|=|Z|=q，。问：M中是否包含一个匹配M，使得（等价问题是求最大三维匹配）。

注：这里给出的是标准形式，一般可不必要求|X|=|Y|=|Z|，表示笛卡尔乘积。三维匹配问题是下一章中二维匹配（2DM）问题的推广，2DM是P问题而3DM是NP完全的。一个匹配是指M的一个子集合{(xi,yi,zi)}，xi∈X，yi∈Y，zi∈Z，且当下标不同时，它们分别取三个集合中的不同元素。3DM可作如下形象的解释：记一单身男人集合为X，一单身女人集合为Y，此外还有一个住房集合Z。其间存在一相容关系（例如有些人之间不愿组成家庭，或不愿住某一住房），这样就给出了一个集合，M是由问题给出的，表示了所有可能组合。所求的匹配即组成的一组家庭（包括住房），其中不能有重婚，也不能让不同的两个家庭住进同一住房。53（3）（划分问题）给定一正整集合A={a1,a2,…,an}，问是否存在A的一个子集A’,使得，即是否可将A中的数分成总和相等的两部分。（4）（顶点复盖问题——VC）给定一个图G=（V，E）及一个正整数K≤|V|，问G中是否有不超过K个顶点的复盖。（一个顶点的子集称为G的一个复盖，若它至少包含G中任一边的两个端点之一）。（5）（哈密顿圈问题——HC）见例8.9。

（6）（团问题）给定图G=（V，E）及一正整数K≤|V|，问是否存在图G中的一个团V’，|V’|≥K。（G的一个顶点的子集V’称为G的一个团，若V’中任意两点间都有E中的边相连）。

图8.5表达了上述六个问题及SAT问题之间的多项式转化关系，即SAT∝3SAT，3SAT∝3DM及3SAT∝VC等等。54例2.（婚姻问题）在遥远的地方有一位酋长，他想把三个女儿嫁出去。假定三个女儿为A、B、C，三位求婚者为X、Y、Z。每位求婚者对A、B、C愿出的财礼数视其对她们的喜欢程度而定,（见下表）：

XYZA3526B271028C147

问酋长应如何嫁女，才能获得最多的财礼（从总体上讲，他的女婿最喜欢他的女儿）。55例2.显然是指派问题的实例，但它也可以看成是两分图赋权匹配问题的实例。用三个点表示酋长的三个女儿，将它们放在一边。再用三个点表示求婚者，将它们放在另一边。在有可能结婚的两人之间画一条边，并在边上写上求婚者对这种结婚愿付出的财礼数，得到右图。右图是一个特殊的图，它的顶点可以分成两个子集，只有分属不同子集的点才可能有边相连（但也可以无边），这样的图称为两分图。定义3.(匹配)

图G的一个匹配是指边集E的一个子集M，M中的任意两条边均不具有公共的顶点。容易看出，酋长要解的问题是在上面的两分图中找出一个具有最大权和的匹配，读者不难由此得到一般两分图最大权匹配问题的数学模型。

以上举的是一个P问题，下面来看一个NP难问题56例3

（装箱问题——Bin—packing）有一批待装箱的物品J={p1,…,pn}，pj的长度为l(pj)。现有一批容量为C的箱子（足够数量），要求找到一种装箱方法，使得所用的箱子数最少。

例3是一个一维的装箱问题。例如，我们有一批具有相同长度的钢材，如果想取出几根已知长度的钢料生产某种设备，当然会希望少用几根原始钢材以减少浪费。此时，我们就遇到了一个一维的Bin—packing问题。当我们想从购买来的三夹板上锯出n块已知长、宽的板来制作家具时，遇到的是二维Bin—packing问题。而当我们真正想把一批已知长、宽、高的物体装入具有相同尺寸的箱子时，又遇到了三维的。下面的定理告诉我们，即使是一维的Bin—packing问题也是NP完全的，故二维和三维的Bin—packing问题更不可能是P问题（它们也是NP完全的）。定理1.

（一维）Bin—packing问题是NP—完全的。证明：易见，划分问题可转化为Bin—packing问题。事实上，取，J={c1,…,cn}可划分为两个相等的集合的充要条件是它们可装入两只容量为C的箱中。

57存在两种不同类型的Bin-packing问题。第一类不允许整理，必须按给定的顺序装箱，称为on-line问题；第二类允许先对物品加以整理，然后再装箱，称为off-line问题。（一）on-line排序问题的近似算法1.NF（NextFit）算法步1将P1装入B1中。步2装Pi,(i=2,…,n)。设Bj是具有最大足码的非空箱，若Pi可继续装入Bj则将其装入，否则将Pi装入Bj+1中，即装入下一空箱中（前面的箱子将不再使用）。记NF算法使用的箱子数为NF（L），则：NF(L)≤2OPT(L)，且2是紧的，即不能被减小。582.FF（FirstFit）算法步1将P1装入B1中。步2装Pi,i=2,…,n。找出最小的j使Pi可装入Bj中，将Pi装入其中，在装箱结束前，每一箱子的剩余空间均有可能被继续利用。定理2

设FF(L)为用FF算法将L中的物品装箱所用的箱数，则有FF(L)≤1.7OPT(L)+1，且1.7是紧的，即不能被减小，（证明有较大难度，从略）。3.BF（BestFit）算法步1装P1装入B1中。步2装Pi,i=2,…,n。在能够装下Pi的箱中找出已装得最多（即剩余空间最小）的一只Bj，将Pi装入Bj。定理3

设BF(L)为用BF算法将L中物品装箱所用的箱数，则有BF（L）≤1.7OPT(L)+1，且1.7是紧的，（证明从略）。59

5.逻辑模型60例1

拟将一批尺寸为1×2×4的的商品装入尺寸为6×6×6的正方体包装箱中，问是否存在一种装法，使装入的该商品正好充满包装箱。解

将正方体剖分成27个2×2×2的小正方体，并按下图所示黑白相间地染色。再将每一2×2×2的小正方体剖分成1×1×1的小正方体。易见，27个2×2×2的正方体中，有14个是黑的，13个是白的（或13黑14白），故经两次剖分，共计有112个1×1×1的黑色小正方体和104个1×1×1的白色小正方体。虽然包装箱的体积恰好是商品体积的27倍，但容易看到，不论将商品放置在何处，它都将占据4个黑色和4个白色的1×1×1小正方体的位置，故商品不可能充满包装箱。61

德国著名的艺术家AlbrechtDürer(1471-1521)于1514年曾铸造了一枚名为“MelencotiaI”的铜币。令人奇怪的是在这枚铜币的画面上充满了数学符号、数字及几何图形。这里，我们仅研究铜币右上角的数字问题

例2.Dürer魔方（或幻方）问题62所谓的魔方是指由1~n2这n2个正整数按一定规则排列成的一个n行n列的正方形。n称为此魔方的阶。Dürer魔方：4阶，每一行之和为34，每一列之和为34，对角线（或反对角线）之和是34，每个小方块中的数字之和是34，四个角上的数字加起来也是34什么是Dürer魔方多么奇妙的魔方！铜币铸造时间：1514年

63

构造魔方是一个古老的数学游戏，起初它还和神灵联系在一起，带有深厚的迷信色彩。传说三千二百多年前（公元前2200年），因治水出名皇帝大禹就构造了三阶魔方（被人们称“洛书”），至今还有人把它当作符咒用于某些迷信活动，大约在十五世纪时，魔方传到了西方，著名的科尼利厄斯·阿格里帕（1486-1535）先后构造出了3~9阶的魔方。64如何构造魔方奇数（不妨n=5）阶的情况Step1:在第一行中间写1Step2:每次向右上方移一格依次填按由小到大排列的下一个数，向上移出界时填下一列最后一行的小方格；向右移出界时填第一列上一行的小方格。若下面想填的格已填过数或已达到魔方的右上角时，改填刚才填的格子正下方的小方格，继续Step2直到填完12345678910111213141516171819202122232425偶数阶的情况

偶数阶的魔方可以利用奇数阶魔方拼接而成，拉尔夫·斯特雷奇给出了一种拼接的方法，这里不作详细介绍65五阶没人知道有多少个！！！三阶1个反射和中心旋转生成8个四阶880个反射和中心旋转生成7040个魔方数量随阶数n增长的速度实在是太惊人了！同阶魔方的个数66允许构成魔方的数取任意实数允许取实数

n阶魔方A、B，任意实数α、βαA+βB是n阶魔方具有指定性质的魔方全体构成一个线性空间问题已发生了实质性变化注：刻画一个线性空间只需指出它的维数并求出此线性空间的一组基底松驰问题的讨论67

1在第一行中共有4种取法，为保持上述性质的成立，第二行中的1还有两种取法。当第二行的1也取定后，第三行与第四行的1就完全定位了，故一共可作出8个不同的最简方阵，称之为基本魔方并记之为Q1，…，Q8

仍以4阶方阵为例。令R为行和，C为列和，D为对角线和，S为小方块和定义0-方：R=C=D=S=0定义1-方：R=C=D=S=4R=C=D=S=1的方阵构成的线性空间具有什么样的性质？类似于构造n维欧氏空间的标准基，利用0和1我们来构造一些R=C=D=S=1的最简单的方阵。68

显然，Dürer空间（简称D空间）中任何一个元素都可以用Q1，Q2，…，Q8来线性表示，但它们能否构成D空间的一组基呢？69容易看出：Q1，…，Q8这8个基本方是线性相关的，即至少存在一个Qj，可以通过其它7个基本方的线性组合得到。这8个基本方的地位是等同的，故可不妨设j=8。下面验证Q1，Q2，…，Q7是否线性相关。

令：，即=70等号两边对应元素相比较，得r1=r2=…=r7=0，所以是线性无关是D空间的最小生成集。

令D

即解方程组：

=

解得

D=研究AlbrechtDürer铸造的铜币712004年浙江大学数学建模竞赛（B题）通讯卫星上的开关设置

地面上存在着n个接收站与n个发送站，而在通讯卫星上则设置了若干种开关模式。开关模式可用矩阵P=(pij)来表示，若卫星可接收发送站i发射的信息并将信息传送回地面的接收站j时，矩阵中的元素pij=1，否则pij=0。通讯卫星上的接收发送任务也可以用一个矩阵T=（tij）来表示，其元素tij为需经通讯卫星传递的由i发点发送到j接收点的信息量的传送时间长度。由于技术上的原因，当发送站i在发送给接收站j信息时，它不能同时发送给别的接收站信息；同样，当接收站j在接收发送站i的信息时，也不能同时接收其他发送站发送的信息。你的任务是：72

设计一组开关模式，k=1,…,r（注：r应当尽可能小），使得对任意给定的任务矩阵T，卫星开关设置均能完成要求的发接收任务。设计一个算法，在发接收任务T给出后，可根据你设计的开关模式（k=1,…,r）求出各开关的使用时间λk，使得在完成预定传送任务的前提下使用各开关模式的总时间最短。同样由于技术上的原因，开关模式的总数r有一个上限。当需要传送的任务数数量较大时，仍无法分派任务。请你想一些办法来解决这一困难，（当然，这时你可能要作出一些牺牲，即传送时间可能会增加一些）。73问题及模型问题的标准形式为：在地面上存在着n个收站与n个发战，而在通讯卫星上则设置了若干种开关模式。开关模式可用矩阵P=(pij)来表示，若卫星可接收发射站i发射的信息并将信息传送回地面的接收站j时，矩阵元素pij=1，否则pij=0。通讯卫星的接发任务也可用一矩阵T=（tij）来表示，其元素tij为需经通讯卫星传递的由i发点发送到j接受点的信息量的传送时间长度。问题要求求r并设计一组开关模式Pk，k=1,…,r及模式Pk的使用时间λk，使得在完成预定传送任务的前提下各开关模式使用的总时间最短，即要求求解下面的问题：74例1

设

这是一个有3个发送站与3个接收站的实例，tij在矩阵中已给出，例如由发站1传送到收站1的通讯量为3单位时间等。分析

容易看出，三个发站需传送的时间分别为6、5、5；而三个收站需接收的时间分别为6、3、7。为完成全部传送任务，通讯卫星总传送时间至少应为7单位时间，即的下界为7。由于技术上的原因，当发站i在发送给收站j信息时，它不能同时发送给别的收站信息；同样，当收站j在接收发站i的信息时，也不能同时接收其他发站发送的信息。这一要求说明，任一开关模式Pk应具有以下性质：（1）Pk的每一行中有且只有一个1，每一列中也有且只有一个1；（2）所有的1均位于不同的行列中。满足（1）、（2）的矩阵被称为置换矩阵，n阶置换矩阵Pk共有n!个，当n较大时，我们不可能在通讯卫星上设置这么多种不同的开关模式。因而，为了设计出切实可行的开关模式，我们还得另想办法。75（设计方法1）注意到Pk每行（或列）元素之和均为1，故不管如何指派开关的使用时间（即不论如何取λk），矩阵均具有某些特殊的性质，例如其行和（及列和）均为同一常数。这样的矩阵构成一个线性空间（参见逻辑模型第一节Dürer魔方），为减少开关模式的种类，可取此空间的一组基底作为开关模式。在使用这种开关模式时，无论T的元素tij怎么取，通讯卫星对每一发（收）点的开通时间总和是恒定的。在这种开关模式下，可按如下方式指派各开关模式的使用时间：步1先将T改变为，满足：（1）≥T（2）记，步2用Pk表示，即将分解为（r为空间的维数）76将T化为的方法一般有无穷多种，如可如下化法：令事实上，，（即通讯卫星传送总时间的下界）。令其中用这种方法化例中的T，得到的任一行（或列）中元素之和均为7。77定义1称行和、列%%和均相等的矩阵为双随机矩阵（Doublystochasticmatrix）定理1

（Birkhoff定理，1944）任一n阶双随机矩阵均可写成至多（n－1）2+1个置换矩阵的非线性组合。的分解可如下进行：步1选取由Pij>0可推出>0的置换矩阵P步2确定

步3取，用－代替步4若=0，停；否则，返回步1。例2.

为方便起见，我们来分解一个元素均为非负整数的3阶双随机矩阵，（由Birkhoff定理，r≤5）78解：取，λ=min{1,3,3}=1分解成，再取

因min{5,5,3}=3，又有，取

79于是又有易得分解结果为：80尚需解决的问题是如何求P，使得Pij>0必有。读者不难发现，此问题可以通过求解一个两分图上的最大流（或最大匹配）来实现，计算量为O(n4)，是多项式时间可解的。具体方法为：作一两分图，若，则作边（i,j），令边容量为1，这样，可作出P的充要条件是该最大流问题的最大流量为n。对例9.33，n=3。由于所有，先取

，－P1为

相应的两分图为：于是又可求得81，相应的两分图为：又可得，…，如此下去，直到作不出P为至，由于的特殊性质及Birkhoff定理，上述分解必能在不超过r=(n－1)2+1步内终止。上述开关设计方法要求在通讯卫星上设置(n－1)2+1种不同的开关模式（即Pk），当n稍大时，(n－1)2+1仍显得太大而使得使用