数学建模

1.Lingo中的一个送分题（8分）

例3：某帆船制造公司要决定下两年八个季度的帆船生产量。八个季度的帆船需求量分别是40条、60条、75条、25条、30条、65条、50条、20条，这些需求必须按时满足，既不能提前也不能延后。该公司每季度的正常生产能力是40条帆船，每条帆船的生产费用为400美圆。如果是加班生产的，则每条生产费用为450美圆。帆船跨季度库存的费用为每条20美圆。初始库存是10条帆船。如何生产？

（注：本题当然可以用动态规划方法求解；接下来你将看到，建立优化模型用Lingo

软件求解比较省事。）

解：八个季度的需求量数组记为xq，则 xq=[40,60,75,25,30,65,50,20]. 类似地，

用数组zc, jb, kc分别表示八个季度的正常生产量、加班生产量、季度末库存量。 目标函数是全部费用之和：min8(400zc(i)450jb(i)20kc(i)).

i1, 约束条件：生产能力 zc(i)40,i1,2,..8.；

数量平衡

kc(1)10zc(1)jb(1)xq(1),kc(i)kc(i1)zc(i)jb(i)xq(i),i2,3,...,8.

以上是模型。怎样用Lingo编程呢？ 把下标的范围当作集合，本题的集合是{1,2,3,4,5,6,7,8}；定义在集合上的一个个数组，都分别称为该集合的属性，本题这个集合有四个属性，分别是xq,zc,jb,kc . 先看本题的Lingo程序，再看注解： model: sets:

jihe/1..8/:xq,zc,jb,kc; endsets data:

xq=40,60,75,25,30,65,50,20; enddata

min=@sum(jihe:400*zc+450*jb+20*kc); @for(jihe:zc<=40);

kc(1)=10+zc(1)+jb(1)-xq(1);

@for(jihe(i)|i#gt#1:kc(i)=kc(i-1)+zc(i)+jb(i)-xq(i)); end

注解：（1）一个完整的Lingo程序必然以“model:”开始，以“end”结束。

（2）集合定义部分（从“sets:”到“endsets”）：内容是定义集合及其属性，命令格式为 “ 集合名/集合的全部元素/:全体属性 ”。本程序中jihe/1..8/:xq,zc,jb,kc; ，这里的“1..8”

等同于“1,2,3,4,5,6,7,8” 。

（3）数据输入部分（从“data:”到“enddata”）：内容是输入已知数据。 （4）其它部分：包括目标函数、约束条件。

本程序中，@sum(jihe:400*zc+450*jb+20*kc)，是求和函数，

@for(jihe:zc<=40)，是循环函数，zc的所有分量都不大于40，

@for(jihe(i)|i#gt#1:kc(i)=kc(i-1)+zc(i)+jb(i)-xq(i))，表示对集合jihe中的所有大于1的i，都满足该项约束。

三．Lingo中的基本集合与派生集合 例4：料场选址问题

六个建筑工地的位置（用平面坐标a、b表示，距离单位：km）及其对水泥的日用量（用d表示，单位：t）由下表给出：

工地的位置（a，b）及 水泥日用量d

-------------------------------------------------------------------------------------------------- 工地编号 1 2 3 4 5 6 a 1.25 8.75 0.5 5.75 3 7.25 b 1.25 0.75 4.75 5 6.5 7.75 d 3 5 4 7 6 11

--------------------------------------------------------------------------------------------------

现有两个临时料场位于P(5,1), Q(2,7)，每日提供水泥的最大能力分别为20t. 假设从料场到各工地均有直线道路连接，运输费用与距离、重量成正比。 （1）请制定运输计划，使总运费尽量低。

（2）进一步调整这两个临时料场的位置，使总运费最低。

解：第i号工地：位置(ai,bi)，水泥日用量di， i=1,2,3,4,5,6.

第j号料场：位置(xj,yj)，水泥日供应能力ej， j=1,2. 从j号料场向i号工地的日运输水泥量记为 cij.

注意：在问题（1）中，(xj,yj)为已知数据，故决策变量为cij，共12个； 在问题（2）中，(xj,yj)待定，故决策变量为xj,yj,cij，共16个。

从j号料场到i号工地，距离为(xjai)2(yjbi)2，送去重量为cij的水泥，二者乘积即为运输费。 目标函数： mincj1i126ij(xjai)2(yjbi)2

2约束条件： 满足需求：

cj1ijdi,i1,2,...,6

供应能力：

ci16ijej，j=1,2

非负性： cij0

这就是本题的优化模型。

尝试用Lingo求解该模型时，六个建筑工地作为一个集合gdjh，两个料场作为一个集合lcjh。接下来就会遇到困难：决策变量cij不仅是依赖于集合gdjh的属性，而且是依赖于集合lcjh的属性，这样的属性应该如何定义呢？

根据两个基本集合gdjh与lcjh构造一个派生集合gdlcjh，再把cij定义为这个集合的属性。

先看本题第(1)问的Lingo程序，再看注解： model: sets:

gdjh/1..6/:a,b,d; lcjh/1,2/:x,y,e; gdlcjh(gdjh,lcjh):c; endsets data:

a=1.25,8.75,0.5,5.75,3,7.25; b=1.25,0.75,4.75,5,6.5,7.75; d=3,5,4,7,6,11;

x,y=5,1,2,7;e=20,20; enddata

min=@sum(gdlcjh(i,j):c(i,j)*((x(j)-a(i))^2+(y(j)-b(i))^2)^0.5); @for(gdjh(i):@sum(lcjh(j):c(i,j))=d(i)); @for(lcjh(j):@sum(gdjh(i):c(i,j))<=e(j)); end

注解：（1）定义派生集合及其属性的命令格式为

派生集合名（基本集合1，基本集合2）：属性

（2）赋值语句“x,y=5,1,2,7”的赋值顺序是“x(1)=5,y(1)=1,x(2)=2,y(2)=7”，而不是“x(1),x(2),y(1),y(2)”；在程序中，该语句可换成“x=5,2;y=1,7”，功能相同。

（3）当表达式中出现的下标符号多于1个时，必须指明针对哪个符号做运算。

2.最小二乘法做数据拟合（10分）

先按照f(x)建立一个函数文件，文件第一句格式为 function f=文件名(参数组, 节点数组)

再用下面命令格式

参数组=lsqcurvefit(‘函数文件名’，参数组初值，节点数组，函数值数组)

例：对函数y=y(x)测量得下面一组数据：

x： 1 2 3 4 5 6 7 8 9 y：4.54, 4.99, 5.35, 5.65, 5.90, 6.10, 6.26, 6.39, 6.50

用函数ya1a2ea3x作拟合，并画图显示拟合效果，给出节点处总误差。

解：先建立并保存函数文件：文件名syp78hswj 内容为：

function f=syp78hswj(a,x0) f=a(1)+a(2)*exp(-a(3)*x0); 再下面主程序： clear hold on

x0=1:9;y0=[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50]; for i=1:9

plot(x0(i),y0(i),'+') end

cscz=[1,1,1];

a=lsqcurvefit('syp78hswj',cscz,x0,y0) x=0:0.2:10;f=syp78hswj(a,x);plot(x,f) f=syp78hswj(a,x0); wc=sqrt(sum((f-y0).^2)) hold off

执行得：(a1,a2,a3)(6.9805,2.9911,0.2031) 节点处总误差 wc=0.0076

问题“录像机中记数器用途”中的参数估计

已知：

t： 0 20 40 60 80 100 120 140 160 183.5 n：0 1153 2045 2800 3466 4068 4621 5135 5619 6152 现在用“最小二乘法”估计参数a,b，使得 tanbn . 先建立函数文件：

function t=lxjjsqhswj(a,n0) t=a(1)*n0.^2+a(2)*n0; 再做主程序： clear

n0=[0,1153,2045,2800,3466,4068,4621,5135,5619,6152]; t0=[0,20,40,60,80,100,120,140,160,183.5]; cscz=[0.01,0.01];

a=lsqcurvefit('lxjjsqhswj',cscz,n0,t0) 执行得：

a = 0.00000250482782 0.01440512449034

23.解微分方程（8分）

数值解

（1）一阶微分方程

2xdyy,0x1,dxy 现以步长h=0.1用“4阶龙格—库塔公式”求数值解： y(0)1.先建立“函数M—文件”：function f=eqs1(x,y)

f=y-2*x/y;

再命令： 格式为：

[自变量,因变量]=ode45(‘函数文件名’,节点数组,初始值) 命令为： [x,y]=ode45('eqs1',0:0.1:1,1) 若还要画图，就继续命令： plot(x,y) （2）一阶微分方程组

cosx2y1y2,0x1,y1sinxy12y2, 只须向量化，即可用前面方法： y2y(0)0.2,y(0)0.3.21function f=eqs2(x,y)

f=[cos(x)+2*y(1)-y(2);sin(x)-y(1)+2*y(2)];

将此函数文件，以文件名eqs2保存后，再下命令： [x,y]=ode45('eqs2',0:0.1:1,[0.2;0.3])

（注：输出的y是矩阵，第i列为函数yi的数值解） 要画图，就命令plot(x,y(:,1)), 及命令plot(x,y(:,2)) （3）高阶微分方程

先化成一阶微分方程组，再用前面方法。

y2yy4x0,0x1,上机练习：

y(0)1,y(0)0.2,y(0)0.准备：令y1y,y2y,y3y，化成

yy21y1(0)1y3y2,y2(0)0.2.

yy2y4xy(0)023134.Monte Carlo法计算面积体积（12分）

已知：yf(x),axb,0f(x)M

求：yf(x),y0,xa,ab所围曲边梯形面积

Monte Carlo法：在矩形axb,0yM内，均匀分布地随机产生n个（大

数）点，落在曲边梯形内的点个数记为m，则

“曲边梯形面积”比“矩形面积”== m比n 例：求ysinx,0x与ox轴所围面积。

分析准备：准确面积==sinxdxcosx02.

0 现用Monte Carlo法计算，取一个矩形为

0x,0y2。

程序：（文件名：syp113） clear

for k=1:1000

a=rand(2,1000);

a(1,:)=a(1,:)*pi;a(2,:)=a(2,:)*2; m=0;

for i=1:1000

f=sin(a(1,i)); if a(2,i)<=f m=m+1; end end

s(k)=2*pi*m/1000; end

ss=mean(s),wc=ss-2

结果： 近似值=？ 误差=？

题：求一个单位球的八分之一部分的体积。

41准确值==r20.52359877559830

386用Monte Carlo法计算：

八分之一单位球：x2y2z21,x,y,z0 一个大正方体：0x,y,z1的体积是1

在“大正方体”内均匀分布产生n个点，记录下落在“八分之一单位球”内的点的个

数m，则m/n即为所求。

程序：（文件名：syp114lx） clear

for k=1:100

a=rand(3,100); m=0;

for i=1:100

f=sum(a(:,i).^2); if f<=1

m=m+1;

end end

s(k)=m/100; end

zqz=pi/6;

ss=mean(s),wc=ss-zqz

结果： 近似值=？ 误差=？

学生练习：

1．求y1exln(2xx2)sinx,0x1与ox轴所围面积。

2．求两个曲面z(1x2y2)sin(exy)与z5x2y2在区域0x,y1上所为立体体积。

二、建模及最优解问题（32分）

1.微分方程模型

例1：某个动物种群，当其数量不多且环境资源丰富时，增长率（即：单位

时间内增长的个体数量与总数量之比）为常数，（这是 马尔萨斯（Malthus）英 1766—1834 提出来的）。为了计算该种群的数量随时间的变化规律而建立的模型称为“马尔萨斯人口模型”。

基本假设：个体数量增长率（即：单位时间内增长的个体数量与总数量之比）为

常数。

模型建立：tt0时刻的数量为已知，记为p0，增长率常数记为k，个体数量p是时间t的函数，pp(t)，在短时间ttt内，增长的数量为

p(tt)p(t)kp(t)t，

得

dp(t)kp(t),p(t0)p0，（这是一个微分方程模型） dt解得 p(t)p0ek(tt0).

例2：一个动物种群，当其数量较多时，增长率显然要受到环境资源的限

制。

基本假设：（这是 逻辑斯谛 提出来的）受资源限制，数量p有最大值M，增长率k

不再是常数，而是随数量p的增大而线性递减，且当p达到最大值M时k=0，记 k=r(M-p)，r为常数 .

模型建立：tt0时刻的数量为已知，记为p0，任意t时刻的数量为p(t). 在短时间ttt内，增长的数量为

p(tt)p(t)kp(t)tr(Mp(t))p(t)t，

得

dp(t)r(Mp)p,p(t0)p0， dt（这是个微分方程模型，称为逻辑斯谛人口模型）

解得 p(t)Mp0 。

p0(Mp0)erM(tt0) 例3：弱肉强食的双种群模型

一个孤立的海岛上：青草；兔子；狐狸。

青草足够多，兔子吃草大量繁殖，兔子多则狐狸容易捕食，狐狸数量增加，导致兔子减少，狐狸也减少。这样，兔子、狐狸数量交替增减，无休止的循环，形成生态平衡。 数学家Volterra（意）建立了这样的“双种群生态模型”：

基本假设：（1）兔子的繁殖率为常数；（2）任何一只兔子与任何一只狐狸相遇的概率是常

数，相遇时兔子被杀的概率是常数；（3）狐狸的死亡率为常数；（4）狐狸的繁殖率与“兔子数乘狐狸数”成正比。

建立模型：t时刻兔子数x(t)，狐狸数y(t)，据基本假设得

dx(t)rx(t)ax(t)y(t),dt dy(t)sy(t)bx(t)y(t).dt 例4：小型试验火箭

初始质量1400kg，其中包括1080kg燃料，竖直向上发射，燃料燃烧速度

18kg/s，产生的推力32000N，空气阻力与速度平方成正比(比例系数0.4kg/m) . 只考虑从点火到燃料燃尽期间的情况，求速度、高度等指标，只要求建立模型。

解： 时间t 高度h(t) 速度v(t)

加速度a(t) 质量m(t) 推力F

空气阻力f(t) 比例系数k 重力加速度g

初始t=0时刻：m(0)=1400, v(0)=0, h(0)=0.

任意t时刻：F=32000,k=0.4,g=9.8,m(t)=1400-18t,

v(t)=h(t), a(t)=h(t), fkv2.

模型的建立：根据“牛顿第二定律”，有

Fmgfma.

模型一：以速度v(t)作为未知函数，建立微分方程模型，

k32000dv2vg,140018t140018t dt0t60.v(0)0,模型二：以高度h(t)作为未知函数，建立微分方程模型，

(140018t)h(t)k(h(t))218gt1400g320000,例5：慢跑者与狗 0t60.h(0)0,h(0)0,一个慢跑者在平面上沿椭圆以恒定的速率v=1跑步,设椭圆方程为

x=10+20cost, y=20+15sint．突然，有一只狗攻击他，他在点（30,20）处，这狗从原点出发,以恒定速率w跑向慢跑者,狗的运动方向始终指向慢跑者．为了分别求出w=20,w=5时狗运动轨迹，请你建立数学模型。

解：

1．人：速度v=1 初始位置（30,20） 逆时针旋转

轨迹椭圆 X(t)=10+20cost, Y(t)=20+15sint

(X10)2(Y20)21 即 222015设t 时刻慢跑者的坐标为(X(t),Y(t))，则方程求导 2(X10)2(Y20)Y0 4002259X10切线斜率 Y

16Y20所以，切线方向为 { -16(Y-20), 9(X-10) }

（思考：若是顺时针，切线方向如何？）

可得 X(t)116(Y20)(16(Y20))(9(X10))9(X10)(16(Y20))(9(X10))2222

Y(t)1

X(0)=30 , Y(0)=20 .

2．狗： 速度w 初始位置（0, 0） 方向指向人

设t 时刻，狗的坐标为(x(t)，y(t))．则 速度方向为

{X(t)-x(t), Y(t)-y(t)}

可得 x(t)wX(t)x(t)(X(t)x(t))(Y(t)y(t))Y(t)y(t)(X(t)x(t))(Y(t)y(t))2222

y(t)w

x(0)=0 , y(0)=0 .

2.画图计算（图论中求定点到其余个点最短路的迪克斯特里算法、动态规划中的图解法 图论中的行遍性最佳邮递员路）

例：求最佳邮递路线，G(V,E),V={v1,v2,,v9},

v1v2(4),v1v7(1),v2v3(2),v2v4(5),v2v7(2),v3v4(1),E={v4v5(1),v4v8(3),v4v9(8),v5v6(10),v6v7(1),v7v8(9),}

v7v9(6),v8v9(3) 解：4个顶点{v4,v7,v8,v9}是奇次，先求他们之间的最短路

P(v4,v7)v4v3v2v7,d(v4,v7)5P(v4,v8)v4v8,P(v4,v9)v4v8v9,P(v7,v8)v7v8,P(v7,v9)v7v9,P(v8,v9)v8v9,d(v4,v8)3d(v4,v9)6d(v7,v8)9d(v7,v9)6d(v8,v9)3

再求这4个点之间的最佳配对，当4与7配对、8与9配对时，距离之和达到最小5+3=8.

在G中添加路径v4v3v2v7与v8v9，形成一个欧拉图，该图的任意一条欧拉路都是最佳邮递路线。

题：求最佳邮递路线

3 2 3 2 2 4 4 3 3 3 3 1 1 3 2 3

三、Matlab中矩阵、数组、rand运算、函数（指数、对数）plot画图（20分）

概念：建模方法、插值、方程求根、图论等（10分）

1. Lagrange插值 定义1

lj(x)(xx0)(xxj1)(xxj1)(xxn)(xjx0)(xjxj1)(xjxj1)(xjxn)

(其中，j分别取0,1,2,,n)称为Lagrange插值基函数，简称基函数。

定义2 称 L(x)yjlj(x) 为Lagrange插值多项式。

j0n注：L(x)就可以作为原插值问题的答案，这种解法称为Lagrange插值法。 用Matlab做Lagrange插值计算：

将定义1、定义2编写成Matlab程序，程序清单如下：

function y=lagr1(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k

p=p*(z-x0(j))/(x0(k)-x0(j)); end end

s=s+y0(k)*p; end y(i)=s; end

该程序文件名记为“lagr1.m”，存放在Matlab工作区。

例 x 0 3 5 7 2.0 9 2.1 11 2.0 12 1.8 13 14 15 y 0 1.2 1.7 1.2 1.0 1.6 求当x分别取0, 0.2, 0.4, 0.6, 0.8, ……, 15时，对应的函数值y，并画出函数y=y(x)的图象.

用Matlab作计算：

x0=[0,3,5,7,9,11,12,13,14,15];

y0=[0,1.2,1.7,2,2.1,2,1.8,1.2,1,1.6]; x=0:0.2:15;

y=lagr1(x0,y0,x) (注释：输出所要求的那76个函数值) plot(x,y) (注释：以x为横坐标、y为纵坐标画出图象) 执行结果：

2. 分段线性插值

3. 三次样条插值

定义3 若有一个分段函数S(x)，满足下面三条：

(1). 在每个小区间 [xi,xi1](i0,1,2,,n1) 上是一个三次多项式； (2). 在每个分段点 xi(i1,2,,n1) 处二阶连续可导； (3). S(xi)yi(i0,1,2,,n) .

x:x0x1xn则称S(x)是插值问题的三次样条插值函数。

y:yyy01n用Matlab做三次样条插值计算： 输入数据：

x0[x0,x1,,xn];y0[y0,];xx0:(小步长):xn;

算函数值：y=spline(x0,y0,x)

画出图象：plot(x,y)

如：上面例2中，用三次样条插值法确定函数y=y(x)，该函数图象为

例3 对函数y1,x[5,5]， 21x令节点xi为 -5、-3、-1、1、3、5，令yi1， 21xi分别用Lagrange法、分段线性法、三次样条法作插值计算，并画图象

解 程序清单为： hold on

x0=[-5,-3,-1,1,3,5];y0=1./(1+x0.^2); x=-5:0.1:5; y1=1./(1+x.^2); y2=lagr1(x0,y0,x); y3=spline(x0,y0,x);

plot(x,y1),plot(x,y2),plot(x0,y0),plot(x,y3) hold off 运行结果：