附录 数学建模实验 数学建模活动是最早出现的数学实验教学模式, 谈到这一话题不得不涉及 "数学建模竞 赛" .中国大学生数学建模竞赛始于一九九二年的北京、西安等八大城市联赛,从一九九四 年起,这一竞赛由教育部高教司和中国工业与应用数学学会主办,每年九月底如期进行.竞 赛的主要目的是模拟大学生毕业后的工作环境, 培养大学生的综合能力 (包括建立数学模型 的能力,计算机应用能力,以及创新能力等等) .在竞赛过程中,参赛同学可以使用计算机 和包括计算机软件在内的一切资料,三人一队密切配合、协同作战,如同参加实际的科研项 目一样解决竞赛问题.在第四天提交一分科研论文(即竞赛答卷) .获奖等级有国家一等奖、 二等奖、省一等奖、省二等奖和成功参赛奖. 围绕大学生数学建模竞赛开展的一系列教学活动, 将数学方法、 知识和一些真实问题联 系起来,将大学生的理想与现实世界和数学联系起来.所谓"数学建模" ,是指人们在解决 实际问题时, 考虑各种逻辑可能性并区分出重要的和次要的因素, 将实际问题抽象成一个数 学模型的过程.对于已经建立的数学模型,可使用计算机求解,得到结果后,要解释、分析、 验证.如果与实际情形相违背,还要修改数学模型重复上述过程.在建立数学模型时,可以 利用因特尔网上丰富的信息资源了解问题的实际背景、 技术资料, 以及科技人员研究此类问 题的现状. 数学建模的模式可以描述为: 面对实际问题→建立数学模型→设计算法→上机计 算→分析结果→用计算机写作科技论文. 正是数学建模活动的开展, 促进了计算机和数学软件在数学教学中的使用, 使得数学实 验课能广泛开展. 这里对数学建模的两个典型问题作了介绍, 并列出五个数学建模实验的练 习题目. 一、蠓虫分类问题 生物学家试图对两类蠓虫(Af 与Apf)进行鉴别,依据的资料是蠓虫的触角和翅膀的长度, 已经测得9只Af和6只Apf的数据(触角长度用x表示,翅膀长度用y表示) Af 数据 x 1.24 1.36 1.38 1.38 1.38 1.40 1.48 1.54 1.56 y 1.27 1.74 1.64 1.82 1.90 1.70 1.82 1.82 2.08 Apf 数据 x 1.14 1.18 1.20 1.26 1.28 1.30 y 1.78 1.96 1.86 2.00 2.00 1.96 现需要解决三个问题: (1)如 何凭借原始资料(Af和Apf的已知数 据被称之为学习样本)制定一种方法, 正确区分两类蠓虫; (2)依据确立 的方法,对题目提供的三个样本: (1.24,1.80),(1.28,1.84),(1.40,2.04) 加以识别; (3)设Af是宝贵的传粉 益虫,Apf是某种疾病的载体, 是否应 该修改分类方法. 1、 问题分析与数学模型 这是一个判别问题, 建模的目标是寻找一种方法对题目提供的三个样本进行判别. 首先 根据学习样本的15对数据画出散点图,图中,Af 用?标记,Apf 用*标记.观察图形, 2.2 93 可以发现,Af 的点集中在图中右下方,而Apf 的点集中在图中左上方.客观上存在一条 直线 L 将两类点分开之间,如果确定了直线L并将它作为 Af 和Apf 分界线,就有了判 别的方法.确定直线 L 应依据问题所给的数据,即学习样本.设直线的方程为 0 3 2 1 = + + w y w x w 对于平面上任意一点 P x y ( , ) ,如果该点在直线上,将其坐标代入直线方程则使方程 成为恒等式,即使方程左端恒为零;如果点 P x y ( , ) 不在直线上,将其坐标代入直线方程, 则方程左端不为零.由于 Af 和Apf 的散点都不在所求的直线上,故将问题所提供的数据 代入直线方程左端应该得到表达式的值大于零或者小于零两种不同的结果. 为了建立判别系统,引入判别函数g(P),当Pxy(,)属于 Af 类时,有g(P)> 0 否则 g(P) < 0 为了对判别系统引入学习机制,在学习过程中将两种不同的状态,以"1"和"-1"表示.当Pxy(,)属于 Af 类时,g(P) = 1,否则 g(P) = —1.取321)(wywxwPg++=(1) 于是由所给数据形成约束条件,这是关于判别函数中的三个待定系数 w1,w2,w3 的线性 方程组: ? ? ? = ? = + + = = + + ) 15 , , 10 ( , 1 ) 9 , , 2 , 1 ( , 1 3 2 1 3 2 1 L L j w y w x w j w y w x w j j j j 这是包括三个未知数w1、w2、w3共15个方程的超定方程组,可以求方程组一种广义解,即 最小二乘解. 2、问题求解与数据结果 根据上面分析写出超定方程组 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = + + ? = + + ? = + + ? = + + ? = + + ? = + + = + + = + + = + + = + + = + + = + + = + + = + + = + + 1 96 . 1 30 . 1 1 00 . 2 28 . 1 1 00 . 2 26 . 1 1 86 . 1 20 . 1 1 96 . 1 18 . 1 1 78 . 1 14 . 1 1 08 . 2 56 . 1 1 82 . 1 54 . 1 1 82 . 1 48 . 1 1 70 . 1 40 . 1 1 90 . 1 38 . 1 1 82 . 1 38 . 1 1 64 . 1 38 . 1 1 74 . 1 36 . 1 1 27 . 1 24 . 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w w (2) 94 下面程序输入有关数据,绘制出散点图,求出超定方程组的最小二乘解,最后在散点图中画 出分类直线 xy=[1.24 1.27;1.36 1.74;1.38 1.64;1.38 1.82;1.38 1.90; 1.40 1.70;1.48 1.82;1.54 1.82;1.56 2.08;1.14 1.78; 1.18 1.96;1.20 1.86;1.26 2.00;1.28 2.00;1.30 1.96]; %学习样本数据 z=[1;1;1;1;1;1;1;1;1;-1;-1;-1;-1;-1;-1]; x=xy(: ,1);y=xy(: ,2);x1=x(1:9);y1=y(1:9);x2=x(10:15);y2=y(10:15); plot(x1,y1,'*',x2,y2,'x'),pause %绘制原始数据散点图 A=[x y ones(x)]; w = A \ z %求方程组的最小二乘解 x=1.10:0.02:1.60; y=(— w(1)*x — w(3) )/w(2); %确定分类直线数据 plot(x1,y1,'x',x2,y2,'o',x,y) %在散点图中画分类直线 程序执行后, 从图形窗口将得到上面的学习样本散点和判别直线图, 从命令窗口得到超 定方程组的最小二乘解 1.1 1.2 1.3 1.4 1.5 1.6 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 w = 7.3101 -3.2041 -3.8236 所以直线方程 0 3 2 1 = + + w y w x w 中的三个待定系数分别为 w1 = 7.3101 w2 = -3.2041 w3 = -3.8236 所以判别直线方程为 7.3101 x — 3.2041 y — 3.8236 = 0 判别函数为 95 g(P) = 7.3101 x — 3.2041 y — 3.8236 (3) 将15个学习样本的所有数据依次代入判别函数g(P),可得 k 1 2 3 4 5 6 7 8 9 g(P) 1.1717 0.5430 1.0096 0.4328 0.1765 0.9635 1.1638 1.6024 0.9156 k 10 11 12 13 14 15 g(P) -1.1934 -1.4778 -1.0111 -1.0211 -0.8749 -0.6006 因为前9个g(P)的值为正数,后6个g(P)的值为负数.根据判别函数g(P)定义知,前9个 学习样本为Af类,后6个学习样本为Apf类.这与学习样本本身是一致的. 题目提供了三个样本供判别,它们的数据列表如下 编号 1 2 3 触角长度x 1.24 1.28 1.40 翅膀长度y 1.80 1.84 2.04 将这三组数据代入判别函数 g(P) = 7.3101 x — 3.2041 y — 3.8236 得计算结果列表 k 1 2 3 g(P) -0.5265 -0.3623 -0.1259 所以,由所给数据用判别函数判别三个新蠓虫的类属,均判为 Apf 类. 二、金属切割问题 某些工业部门 (如贵重石材加工等) 采用截断切割的加工方式从一个长方体中加工出一 个已知尺寸、位置预定的长方体(这两个长方体的对应表面是平行的) ,通常要经过六次截 断切割. 已知待加工长方体和成品长方体的长、宽、高分别为 10、14.5、19 和3、2、4, 二 者左侧面、正面、底面之间的距离分别为 6、7、9 (单位均为厘米). 切割费用为每平方厘米 1 元. 试为这些部门设计一种安排各面加工次序的方案,使加工费用最少. 1. 问题分析与数学模型 首先考虑各种不同切割方式的数学描述. 将左、右、前、后、上、下相应的切割面编号 为1、2、3、4、5、6. 则,一个切割方式就是各加工面 { 1, 2, 3, 4, 5, 6 }的一个全排列,记为:σ σ σ σ σ σ 1 2 3 4 5 6 .将所有不同的切割方式组成的集合记为 Σ σ σ σ σ σ σ σ σ σ 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 i i i i j j 现在考虑任一种切割方式σ σ σ σ σ σ σ = ( 1 2 3 4 5 6 ) 的费用的描述. 由于这种切割方 式在操作过程中总共要进行六次切割,记六次切割 的面积依次为 S( ) σ 1 , S( ) σ σ 1 2 , S( ) σ σ σ 1 2 3 , S( ) σ σ σ σ 1 2 3 4 , S( ) σ σ σ σ σ 1 2 3 4 5 ,……………, S( ) σ σ σ σ σ σ 1 2 3 4 5 6 则,六次切割的总费用为 96 f S i i ( ) ( ) σ σ = = ∑ 1 1 6 Lσ (4) 由上此可知,截断切割问题实际上是在可行域 Σ 上求目标函数 ) (σ f 的最小值.由于 这一问题规模不是很大,可以用穷举法求解. 穷举算法 第一步:列出所有可能的切割方案,即720 个操作数(记为一个 720*6 矩阵); 第二步:计算每一种切割方案的总费用(记为 720 个元素的向量); 第三步:从720 个费用数据中选出最小值; 第四步:列出费用最少的全部操作方案,结束. 在这一算法中,比较困难的是第二步的算法实现. 现将有关细节考虑如下: (1)待加工长方体的长、宽、高为 10、14.5、19 用向量记为 l0 10 14 5 19 = ( . ) ) (2)成品长方体的长、宽、高为 3、2、4,它距待加工长方体左侧面、正面、底面之 间的距离分别为 6、7、9. 所以,切割时左,右,前,后,上,下各面的切割厚度数据用向 量记为 h = ( . 6 1 7 55 6 9 (3)h 有六个数据,在切割过程中,将待加工长方体切割为成品长方体要经历六种状 态,每一状态是从前一状态长方体的长、宽、高数据中的某一个减去对应的切割厚度数据而 形成新的长方体. 被减数只有三个而减数有六个,其.对应关系列表为 对应 对应 对应 被减数指标 1 2 3 减数指标 1,2 3,4 5,6 例如当切割左面或切割右面时,都应该从长、宽、高数据中的第一个数减去切割厚度数据. (4) 切割面积的计算. 每一次切割面的面积实际上是本次切割状态中半成品长、宽、 高数据中不改变的两个数据的乘积. 穷举法程序 l1=[1 2 3 4 5 6];k=0; %设置 l1 1 2 3 4 5 6 for j1=1:6,v1=l1(j1);l2=l1;l2(j1)=[]; %从 中选取 v1 并设置 l l1 2 for j2=1:5,v2=l2(j2);l3=l2;l3(j2)=[]; %从 中选取 v2 并设置 l2 l3 for j3=1:4,v3=l3(j3);l4=l3;l4(j3)=[]; %从 中选取 v3 并设置 l l3 4 for j4=1:3,v4=l4(j4);l5=l4;l5(j4)=[]; %从 中选取 v4 并设置 l4 l5 for j5=1:2,v5=l5(j5);v6=l5(3-j5); %从 中选取 v5 与v6 l5 k=k+1;p(k,:)=[v1 v2 v3 v4 v5 v6]; %将[v1 v2 v3 v4 v5 v6]赋于 P end,end,end,end,end 97 clear l1 l2 l3 l4 l5 j1 j2 j3 j4 j5 v1 v2 v3 v4 v5 v6 l0=[10 14.5 19];h=[6 1 7 5.5 6 9];u=[1 1 2 2 3 3]; for k=1:720 v=p(k,:);l=l0;s=0; %提取 P 的第 k 行数据 for j=1:6 i=v(j);l(u(i))=l(u(i))-h(i); %模拟截割第 i 面并减去相应的厚度 ll=l;ll(u(i))=[];s=s+ll(1)*ll(2); %累加计算本次截割面的费用 end q(k)=s; %将六次截割的总面积赋值给 Q end clear h j i k l l0 ll s u v R an=min(q);find(q==an);p(ans,:), %寻求 Q 中的最小者 an 运行程序 3,操作数有两组 6 3 1 5 4 2 6 3 5 1 4 2 使得切割费用为: 374 达到最小值. 对应的切割方案列表如下 费用最少的切割方案 所需切割费用 底面,前面,左面,上面,后面,右面 底面,前面,上面,左面,后面,右面 374 374 三、建模实验练习题 1、 拐角问题 考虑下面问题,分别建立数学模型并求解.讨论三个问题的差异,是否有相同之处. (1)紧靠楼房的温室伸入花园 a=2 米,高b=3 米,温室正上方是楼房的窗台.清洁工 打扫窗台周围时, 需要将梯子从花园地上越 过温室上方靠在楼房墙上. 满足要求的梯子 最小长度应该是多少? x α x α (2)在一个煤矿的矿道里,矿工扛了一架梯 子打算通过交叉矿道口.已知横向矿道宽 2 米,纵向矿道宽 1.5 米.矿道相交成角度 120o .问通过矿道的梯子限长多少? (3)在医院走廊上, 一个立方体的医疗设备要水平运送到手术室,途经直角拐角处.如果横向 走廊宽 2 米,纵向走廊宽 1.2 米.问该设备的长和宽的限定长度为多少? 手术室 手术室 α α 98 2、 捕食者与被捕食者(狐狸与野兔)模型 假设狐狸和野兔共同生活在同一个有限区域内, 有足够多的食物供野兔享用, 而狐狸仅以野 兔为食物.x 为野兔量,y 表狐狸数量.假定在没有狐狸的情况下,野兔增长率为 400%. 如果没有野兔,狐狸将被饿死,死亡率为 90%.狐狸与兔子相互作用的关系是,狐狸的存 在使兔子受到威胁,且狐狸越多兔子增长受到阻碍越大,设增长的减小与狐狸总数成正比, 比例系数为 0.02. 而兔子的存在又为狐狸提供食物, 设狐狸的单位时间的死亡率的减少与兔 子的数量成正比,设比例系数为 0.001.因此数学模型为 ? ? ? ? = ′ + ? = ′ xy x x xy y y 02 . 0 4 001 . 0 9 . 0 (1) 以x=800,y = 100 为初始值,计算 x(t),y(t)当t∈[0,14]时的数据.绘出下面图 形并分析捕食者和被捕食者的数量变化规律. 0 2 4 6 8 10 12 14 0 5 00 10 00 30 00 (2)以x为横坐标,y 为纵坐标绘制如下相位图 15 00 20 00 25 00 兔子数量数量狐狸数量时间99 根据图形分析被捕食者数量增加(或减少)对捕食者数量的影响. 3、 抛物体运动的数学模型 一般的抛物体运动的数学模型为 ? ? ? ? ? ? = = 2 2 1 sin cos gt t v y t v x α α 其中,v 是初速度,α 是发射角,t 是时间参数.考虑如下问题: (1)一般的抛射线运动中的远程炮最大射程; (2)消防水龙头喷水灭三层楼房窗口内火焰时,高压龙头喷射水柱的最远距离; (3)投蓝球的最佳角度; (4)掷铅球的最佳角度; (5)对几个问题的抽象——一般的数学模型 4、 海浪高度数据的拟合问题 海洋水文观测站记录了某海域每天 24 小时海浪潮高度数据(相对于海堤上的零标尺记 号) .工作人员每小时做一次记录,下面是从 12 月1日早上零点开始到 12 月2日晚上 23 点结束两天内完整的数据记录.能否根据这些数据预测三天以后(即12 月5日)的海浪高 度数据 12 月1日海浪高度数据(单位:米) 0:00 1:00 2:00 3:00 4:00 5:00 6:00 7:00 2.4000 1.2000 -0.1000 -0.5000 -2.5000 -3.0000 -2.7000 -1.6000 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 0.2000 2.1000 3.4000 3.6000 2.9000 1.6000 0.2000 -1.2000 16:00 17:00 18:00 19:00 20:00 21:00 22:00 23:00 -2.4000 -3.0000 -3.1000 -2.3000 -0.7000 1.3000 2.9000 3.6000 12 月2日海浪高度数据(单位:米) 0:00 1:00 2:00 3:00 4:00 5:00 6:00 7:00 3.1000 2.0000 0.6000 0.6000 -2.2000 -3.6000 -3.2000 -2.5000 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 -0.9000 -1.1000 2.9000 3.9000 3.6000 2.5000 1.0000 16:00 17:00 18:00 19:00 20:00 21:00 22:00 -2.4000 -3.0000 -3.4000 -3.0000 -1.7000 0.2000 2.2000 3.5000 0 5 0 0 1 0 0 0 1 5 0 0 2 0 0 0 2 5 0 0 3 0 0 0 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 3 5 0 4 0 0 (9 0 0 ,2 0 0 ) 1 2 3 4 狐狸数量兔子数量100 (1) 由数据画散点图 0 1 0 2 0 3 0 4 0 5 0 - 4 - 3 - 2 - 1 0 1 2 3 4 (2) 用快速傅里叶变换(FFT)分析数据,计算出数据变化的频率(或周期)b (3) 以x表示海浪高度变量,考虑用三角函数作数据拟合,设)2cos( ) 2 sin( ) ( 2 1 bt a bt a t x π π + = 利用数据拟合的最小二乘法确定系数a1和a2. (4) 根据确定的数学表达式计算出三天以后(即12 月5日)的海浪高度数据为预测值. (5) 分析拟合误差,修改拟合函数使之更符合实际数据. 5、 渡口渡轮的平均载量问题 渡口的渡船甲板长 32 米,可以并排停放两列车辆.怎样 在甲板上安排过河车辆的位置,才能安全地运过最多数量 的车辆,一次可以运多少辆车. 观察发现每次情况不尽相 同,有下列数据和情况: 车辆随机到达,形成一个等待上船的车列; 车辆中轿车 40%,卡车 55%,摩托 5%; 轿车长为 3.5—5.5 米, 卡车车身长为 8 一10 米, 摩托数量少可安插在车队中间设长度为零. 完成如下练习: 用计算机产生两个随机数R1和R2.第一个随机数R1模拟下一辆到来的是卡车(k=1),轿车(k=2) 或者摩托(k=3).第二个随机数randn模拟到来车辆的长度:卡车L1=9+randn/3, 轿车 L2=4.5+rand/3,摩托L3=0 40% 55% 5% k=1 k=2 k=3 (1)模拟 20 辆车辆到来的情况,列表记录每一辆车到达时的五个项目: 随机数R1;车辆种类;随机数;车身长;队列长 (2) 统计并计算每一次渡轮甲板上能容纳的平均车辆数(卡车、轿车、摩托车分类计算) . 6、 分形曲线的算法实验 Koch 分形曲线由两点 P1、P2 连成的线段开始,将线段分为三段,用一等边三角形去掉 底边后的图形取代中间一段.形成右图所示的图形,这就是 Koch 分形曲线的生成元.由两 点连线演化为五点连线(右图)的算法如下 第一步:Q1 = P1 + (P2-P1)/3; 第二步:Q3 = P1 + 2(P2-P1)/3; 第三步:Q2 = Q1 + (Q3-Q1)*A' ; 101 第四步:P5 ← P2;P2 ← Q1;P3 ← Q2;P4 ← Q3. 第五步:画P1,P2,P3,P4,P5 连线的图形,结束. 在算法中,符号"←"表示右边数据赋值给左边的变量,P1表示第一个点的坐标(x1,y1) , P2表示第二个点坐标(x2,y2)其余类推. 参考下面程序 p=[0 0;10 0]; a=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)]; for k=1:5 n=max(size(p));d=diff(p)/3; q=p(1:n-1,:);p(5:4:4*n-3,:)=p(2:n,:); p(2:4:4*n-6,:)=q+d; p(3:4:4*n-5,:)=p(2:4:4*n-6,:)+d*a'; p(4:4:4*n-4,:)=q+2*d; 分形第二种生成元 end plot(p(:,1),p(:,2)) 除了 Koch 分形曲线的生成元, 还有右图所示的第二 种生成元和第三种生成元.试分析三种生成元的构 成,设计算法,完成以下作业 (1)从一个任意三角形(或四边形)出发,绘制封 闭的 koch 分形曲线; (2)以第二种生成元绘制分形曲线; 分形的第三种生成元 (3)以第三种生成元绘制分形曲线; (4) 修改第三种生成元为树枝形态, 绘树状分形曲 线. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1 0.9 102