数学模型

学生实验报告

实习时间：2015-- 2016学年第二学期专业班级：信息与计算科学1402班

姓名（学号）：李潇亚（20145358）

2016年 5月 8 日

实验名称实验一：人口预测和控制

实验地点实验日期

学时 2

一、实验目的

1.了解最小二乘拟合的基本原理和方法；

2.掌握matlab作最小二乘多项式拟合和曲线拟合的方法

3.通过实例学习如何用拟合方法解决实际问题。

4.了解各种参数辨识的原理和方法；

5.通过范例展现由机理分析确定模型结构，拟合方法辨识参数，误差分析等求解

实际问题的过程；

通过该实验的学习，掌握几种基本的参数辨识方法，了解拟合的几种典型

应用，观察不同方法得出的模型的准确程度，学习参数的误差分析，进一步了解

数学建模过程。这对于学习深入理解数学概念，掌握数学的思维方法，熟悉处理

大量的工程计算问题的方法具有十分重要的意义。

二、实验题目

对美国1790-2000年的人口数据选择合理的方法进行拟合和预测，具体要求如下：

模型1

(1) 假设美国人口上限为4亿，假设1790年的人口增长率为2.95%，计算每隔

十年的人口增长率，并进行适当的处理，建立微分方程模型；

(2) 利用 (1) 中的模型计算各年人口，与实际人口数量比较，计算模型的计算

误差；

(3) 利用 (1) 中的模型预测美国2010,2020,2030,2040,2050年的人口；

模型2.

（1) 根据表中的人口数据，进行曲线拟合，建立数学模型；

(2) 利用 (1) 中的模型计算各年人口，与实际人口数量比较，计算模型的计算

误差；

(3) 利用 (1) 中的模型预测美国2010,2020,2030,2040,2050年的人口；

数据表1 美国人口

年1790 1800 1810 1820 1830 1840 1850 1860 人口/百 3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4

万

年1870 1880 1890 1900 1910 1920 1930 1940

38.6 50.2 62.9 76 92 106.5 123.2 131.7 人口/百

万

年1950 1960 1970 1980 1990 2000

150.7 179.3 204 226.5 251.4 281.4

人口/百

万

三、实验过程

源程序：t=1790:10:2000;

x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 ...

92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4 281.4];

%plot(t,x,'*');

x1=log(x);

a=polyfit(t,x1,1)

x2=exp(a(2)+a(1).*t);

t1=2010:10:2050;

plot(t,x2,'+',t,x,'*');

x3=exp(a(2)+a(1).*t1);

er=abs(x2-x)./x

t2=2010;

x4=exp(a(2)+a(1).*t2)

四、实验结果

a =

0.0202 -34.3934

er =

0.2645

x4 =

516.7091

五、结果分析

在MATLAB命令窗口运行该程序，输出结果a = 0.0202 -34.3934；x1 =516.7091因此，参数a=0.0202, b=-34.3934，拟合函数在2010处的函数值f(2010)=516.7091。误差er=0.2645。

由上图1750-2000的实际人口和预测拟合的值发现，在接近2010年时存在较大误差，故预测的值存在不准确性，还需进一步计算，改善程序，得出较为准确的预测值。

实验名称实验二：用MATLAB求解线性规划问题

实验地点实验日期

学时 2

一、实验目的

1．了解线性规划的基本内容

2.熟悉MATLAB软件求解线性规划问题的基本命令

3.学习灵敏分析问题的思维方法

二、实验题目

2.1.奶制品的生产和销售，用MATLAB求解

课后题1

三、实验过程

2.1奶制品的生产和销售

源程序： c=[-72 -64];

A=[1 1;12 8;3 0];

b=[50;480;100];

lb=zeros(2,1);

[x,fval]=linprog(c,A,b,[],[],lb)

Optimization terminated.

课后题 1：（1）根据题意，假设市政，代办机构，政府，政府，市政分别为x1,x2,x3,x4,x5,总利润为z,得：

Max z=4.3%x1+2.7%x2+2.5%x3+2.2%x4+4.5%x5

x2+x3+x4>=400

2*x1+2*x2+x3+x4+5*x5<=7

9*x1+15*x2+4*x3+3*x4+2*x5<=25

X1+x2+x3+x4+x5<=1000

X1,x2,x3,x4,x5>=0

源程序： c=[-0.043 -0.027 -0.025 -0.022 -0.045];

A=[-1 -1 -1 -1 -1;2 2 1 1 5;9 15 4 3 2;1 1 1 1 1];

b=[-400;7;25;1000];

lb=zeros(5,1);

[x,fval]=linprog(c,A,b,[],[],lb)

（3）由第一问的得，c=[-0.045 -0.027 -0.025 -0.022 -0.045];

A=[-1 -1 -1 -1 -1;2 2 1 1 5;9 15 4 3 2;1 1 1 1 1];

b=[-400;7;25;1000];

lb=zeros(5,1);

[x,fval]=linprog(c,A,b,[],[],lb)

c=[-0.043 -0.027 -0.024 -0.022 -0.045];

A=[-1 -1 -1 -1 -1;2 2 1 1 5;9 15 4 3 2;1 1 1 1 1];

b=[-400;7;25;1000];

lb=zeros(5,1);

[x,fval]=linprog(c,A,b,[],[],lb)

四、实验结果

2.1

x =

20.0000

30.0000

fval =

-3.3600e+003

课后题1

（1）x =

0.0000

0.0000

0.0000

527.1516

2.9486

fval =

-11.7300

所以若该经理有1000万元资金，取x4=527.1516,x5=2.9486,所以投资 D政府527.1516万元和E市政2.9486万元。

（3）x =

0.0000

0.0000

0.0000

527.1516

2.9486

fval =

-11.7300

若证券A的税前收益增加为4.5%，投资不改变。

x =

0.0000

0.0000

0.0000

527.1516

2.9486

fval =

-11.7300

若证券C的税前收益减少为4.8%，投资不应该改变。

五、结果分析

2.1奶制品的生产和销售

上面的结果得到了问题的最优解和最优值，即用20桶牛奶生产A1，30桶牛奶生产A2，但是如果目标函数的系数改变，结果也会有改变，而且上面的数值也是近似的满足条件，所以问题应该进一步深入探讨。

课后题1

（1）通过计算的结果x4=527.1516,x5=2.9486，所以应该投资D政府527.1516万元和E市政2.9486万元。但是共有1000万元资金，应该有更好的方案去投资，获得最大的利润。

（2）题目中提到以2.75%的利率街道不超过100万元资金，但是不知道该问跟第一问的1000万元是否还有关系，所以没有做出来。

（3）分别改变两种证券的税前收益，但是结果不变。

实验名称实验三：用MATLAB求解整数规划问题

实验地点实验日期

学时 2

一、实验目的

1．了解整数规划的基本内容

2.熟悉MATLAB软件求解整数规划问题的基本命令

3.学习灵敏分析问题的思维方法

二、实验题目

课后题2或者3

三、实验过程

课后题2

每个销售代理点只能向本区和一个邻区的大学生售书，可以用0-1变量表示是否建在该区，从而建立这个问题的0-1规划模型。蒋大学生数量为24、29、42、21、56、18、71的区标为1、2、3、4、5、6、7，用rij为第i区的大学生人数，用0-1变量Xij=1表示（i，j）区的大学生由一个代售点供应图书（i

Max=63*x12+76*x13+71*x23+50*x24+85*x25+63*x34+77*x45+39*x46+92*x47+50 *x24+74*x56+89*x67

X12+x13+x23+x24+x25+x34+x45+x46+x47+x56+x67<2

X12+x13<=1

X12+x23+x24+x25<=1

X13+x23+x34<=1

X25+x45+x56<=1

X46+x56+x67<=1

X47+x67<=1

Xij=0或xij=1

源程序：

c=[-63 -76 -71 -50 -85 -63 -77 -39 -92 -74 -89];

>> A=[1 1 1 1 1 1 1 1 1 1 1

1 1 0 0 0 0 0 0 0 0 0

1 0 1 1 1 0 0 0 0 0 0

0 1 1 0 0 1 0 0 0 0 0

0 0 0 1 0 1 1 1 1 0 0

0 0 0 0 1 0 1 0 0 1 0

0 0 0 0 0 0 0 1 0 1 1

0 0 0 0 0 0 0 0 1 0 1];

>> b=[2;1;1;1;1;1;1;1];

>> [x,fval]=bintprog(c,A,b)

四、实验结果

x =

1

1

fval =

-177

五、结果分析

计算得知x25=1,x47=1,所以应该建在2和4区，分别向5区和7区销售，能使所能供应的大学生数量最大，最大为177

实验名称实验四：用MATLAB求解微分方程，传染病模型实验地点实验日期

学时 2

一、实验目的

1.练习数值微分的计算

2.掌握用MATLAB软件求微分方程初值问题的方法

3.通过实例学习用微分方程模型解决简化的实际问题

4.了解欧拉方法和龙格-库塔方法的基本思想和计算公式等概念

二、实验题目

1．指数传播模型

2.SI传播模型

3.SIS传播模型

画出每类人群的图像，相轨线，并讨论阈值的取值对传染病传播速率的影响

三、实验过程

1.指数传播模型

设时刻t的病人人数x(t)是连续、可微函数，并且每天每个病人有效接触（足

以使人致病的接触）的人数为λ，则t到t+Δt病人人数的增加，就有 x(t+Δt)-x(t)=λx(t)Δt

再设t=0时有x0个病人，得：dx/dt=λx,x(0)=x0

2.SI传播模型

设总人数N不变，人群分为易感染者和易感染者两类，记作s(t),i(t),每天每个病人有效接触的人数为λ,为日增长率，初始时刻病人得比例为i0,则

Ndi/dt=λNsi

s(t)+i(t)=1

di/dt=λi(1-i),i(0)=i0

解得：i(t)=1/(1+(1/i0-1)*exp(-λt))

源程序：clc

x0=0.2;

lamta=0.006;

t=0:5:1000;

n=length(t);

x=1+(1/x0-1)*exp(-lamta*t);

x=1./x;

subplot(1,2,1);

plot(t,x,'.-');

legend('x~t曲线');

dx=lamta.*x.*(1-x);

subplot(1,2,2);

plot(x,dx);

legend('dx~x曲线');

3.SIS模型

设每天被治愈的病人数占病人总数的比例为日治愈率常数μ，接触数σ是整个传染期内每个病人有效接触的平均人数，得

Ndi/dt=λNsi-μNi

di/dt=λi(1-i)-μi,i(0)=i0

=λ/μ

解得：di/dt=-λi[i-(1/σ)]

源程序：clear clc

x0=0.1;

lamta=0.1;

sgm=3;

t=0:5:200;

a=1-1/sgm;

x=a+(1/x0-a)*exp(-lamta*a*t);

x=1./x;

subplot(1,2,1);

plot(t,x,'.-');

legend('x~t曲线');

dx=-lamta.*x.*(x-1+1/sgm);

subplot(1,2,2);

plot(x,dx);

legend('dx~t曲线')

四、实验结果

1.结果表明，随着t的增加，病人人数x(t)无限增长，不符合实际

2.SI的i(t)~t和di/dt的图形如图所示，并且由图可知当i=0.5时，di/dt达到最大值，这时病人增加的最快，当t->∞,i->1时，所有人都会被感染，不符合实际。

3.SIS模型中σ是一个阈值，如图所示，当σ>1时i(t)的增加性取决于i0的大小，但其极限值i(∞)=1-1/σ随σ的增加而增加，当σ<=1时，病人比例i(t)越来越小，最终趋于0.

五、结果分析

1.指数模型之所以不成立是因为病人有效接触的人群中，有健康人还有病人，但是只有健康人才会被传染成病人，所以在改进的模型中必须区别这两种人。

2.SI 模型考虑了健康者和病人，得出的结果更加准确了，在计算出来的i=0.5的时候应该加强卫生医疗推迟传染病高潮，但是当当t->∞,i->1时，所有人都会被感染，不符合实际，原因是没有考虑到会有被治愈的人，所以应该更改模型。

3.SIS 模型综合了各种假设，得出最终的结果，但是需要考虑σ的值，当σ>1时i(t)的增加性取决于i0的大小，但其极限值i(∞)=1-1/σ随σ的增加而增加，当σ<=1时，病人比例i(t)越来越小，最终趋于0.

实验名称 实验五：种群问题

实验地点 实验日期

学时 2

一、实验目的

1.练习数值微分的计算

2.掌握用MATLAB 软件求微分方程初值问题的方法

3.通过实例学习用微分方程模型解决简化的实际问题

4.了解欧拉方法和龙格-库塔方法的基本思想和计算公式，平衡点及其稳定性等概念

二、实验题目

两种群相互竞争模型如下：

()1(11)12

()2(12)

12x y x t r x s n n x y y t r y s n n ?

=--???

?=--??

其中x （t ），y(t)分别是甲乙两种群`的数量，r1，r2为它们的固有增长率，n1，n2为它们的最大容量。S1的含义是，对于供养甲的资源而言，单位数量乙（相对n2）的消耗量为单位数量甲（相对n1）消耗的s1倍，对于s2也可做相应的解释。

分析：

这里用x (t)表示甲种群在时刻t的数量，即一定区域内的数量。Y(t)表示乙种群在时刻t的数量。假设甲种群独立生活时的增长率（固有增长率）为r1，则x (t)/ x=r1，而种群乙的存在会使甲的增长率减小，且甲种群数量的增长也会抑制本身数量的增长，即存在种间竞争。这里，我们设增长率的一部分减少量和种群乙的数量与最大容纳量的比值成正比，与s1（s1表示最大容纳量乙消耗的供养甲的资源是最大容纳量甲消耗该资源的s1倍）成正比。另一部分的减少量和种群甲的数量与甲的最大容纳量的比值成正比。则我们可以得到如下模型：

x(t)=r1*x*(1-x/n1-s1*y/n2)

同样，我们可以得到乙种群在t时刻的数量表达式：

y(t)=r2*y*(1-s2*x/n1-y/n2)

如果给定甲、乙种群的初始值，我们就可以知道甲、乙种群数量随时间的演变过程。

问题一：

设r1=r2=1,n1=n1=100,s1=0.5,s2=2, 初值x0=y0=10，计算x(t),y(t)，画出它们的图形及相图(x,y)，说明时间t充分大以后x(t),y(t)的变化趋势（人民今天看到的已经是自然界长期演变的结局）。

问题二：

改变r1,r2,n1,n2,x0,y0，但s1,s2不变，（或保持s1<1,s2>1），计算并分析所得结果；若s1=1.5(>1)，s2=0.7(<1)再分析结果，由此你的得到什么结论，请用各参数生态学上的含义作出解释。

问题三：

实验当s1=0.8(<1),s2=0.7(<1)时会有什么样的结果：当s1=1.5(>1),s2=1.7(>1)时又会有什么样的结果。能解释这些结果吗？

三、实验过程

问题一：

源程序

fun.m:

function dx=fun(t,x)

r1=1;r2=1;n1=100;n2=100;s1=0.5;s2=2

dx=diag([r1*(1-x(1)/n1-s1*x(2)/n2);r2*(1-s2*x(1)/n1-x(2)/n2)])*x;

ts=0:0.1:10;

x0=[10,10];

[t,x]=ode45(@fun,ts,x0);

plot(t,x),grid,

四、实验结果

问题一：

通过上面程序画出三个图关于x1(t),x2(t)的图形，由图可知，两个种群的单个成员消耗的资源差不多相同，而环境能承受的种不同，两个种群竞争结果取决于s1,s2.

五、结果分析

问题一：当s1<1,s2>1时，s1意味着在对供养甲的资源的竞争中乙弱于甲，s2意味着在供养乙的资源的竞争中甲强于乙，于是种群乙终将灭绝，种群甲趋向最大容量。

问题二：当s1>1,s2<1时，s1意味着s1意味着在对供养甲的资源的竞争中乙强于甲，s2意味着在供养乙的资源的竞争中甲弱于乙，于是种群甲终将灭绝，种群甲趋向最大容量。

问题三：当s1<1,s2<1时,因为在竞争甲的资源乙较弱，而在竞争乙的资源中甲较弱，于是可以达到一个双方共存的稳定的平衡状态。

实验名称实验六:综合评价问题

实验地点实验日期

学时 2

一、实验目的

1.掌握各种类型数据的处理方法及MATLAB操作

2.掌握用MATLAB软件求理想解法的基本命令

3.通过实例学习层次分析中一致性的判定及MATLAB求解的基本命令

二、实验题目

校长奖学金评选问题

目前大多数高校都设立了校长奖学金，每年都从各方面表现优秀的在校大学生中选部分最优者作为校长奖学金的获得者。现已知某高校经过层层选拔，推荐除了5名候选者，评委会需要根据各项考核指标从这些学生中选出2名优秀的同学作为获奖对象。考核指标有学分绩点，团队协作能力，体育水平，创新实践能力，现场答辩表现。现在评委会根据推荐单位及候选学生提交的材料，并通过答辩得到了候选者的量化评分结果，见下表：

学分绩点团队协作能

力体育水平创新实践能

力

现场答辩表

现

学生A 95 4 85 90 98 学生B 92 2.5 90 92 95 学生C 94 5 78 88 85 学生D 91 3 83 94 88 学生E 93 3.5 95 92 90

三、实验过程

优秀学生

学分绩点团队协

作能力

体育水

平

创新实践

能力

现场答辩能

力

设学生为i,考核标准为j ，每项用Xij(i<=5,j<=5)表示，假设他们的对比矩阵为A=[1 1/2 1/4 1/2 1/2;2 1 1/2 1 1;4 2 1 2 2;2 1 1/2 1 1;2 2 1/2 1 1]; 源程序：

A=[1 1/2 1/4 1/2 1/2;2 1 1/2 1 1;4 2 1 2 2;2 1 1/2 1 1;2 2 1/2 1 1]; [n,n]=size(A);

RI=[0 0 0.58 0.9 1.12 1.24 1.32 1.41 1.45 1.49 1.51]; [V,D]=eig(A); lam=D(1,1);

CI=(lam-n)/(n-1); CR=CI/RI(n)

计算特征值和特征向量：

A=[1 1/2 1/4 1/2 1/2;2 1 1/2 1 1;4 2 1 2 2;2 1 1/2 1 1;2 2 1/2 1 1]; [x,lumda]=eig(A); r=abs(sum(lumda)); n=find(r==max(r)); max_lumda=lumda(n,n) max_x=x(:,n)

最后， 通过公式1*n

i Bi Aij wi ==∑计算各高校 的最终得分,运用matlab 计算。

四、实验结果

max_lumda =

5.1926

max_x =

-0.1805 -0.3610 -0.7221 -0.3610 -0.4306 ans =

-0.9373

-1.8746

-3.7493

-1.8746

-2.2356

五、结果分析

通过层次分析法，算出这个矩阵通过组合一致性比率，并且通过matlab算出λλ=5.19特征向量wi=[-0.1805;-0.3610;-0.7221;-0.3610;-0.4305].因此最终结果为选择A同学，由于B和C同学得分一样，还需继续计算，得出最终结果。

实验名称实验七:综合实验

实验地点实验日期

学时 4

一、实验目的

1．掌握分析问题和解决问题的过程

2.掌握前几章学习过的利用MATLAB求解方法

3.综合利用所学方法建立模型，并用MATLAB进行求解

二、实验题目

降落伞的选择

为向灾区投放救灾物资2000kg，需选购一些降落伞，已知空投高度为500m，要求降落伞落地的速度不能超过20m/s。降落伞可视为半径为r的半球面，用每根长为l的16根绳索连接着载重m，如图1所示。

每个降落伞的价格由三部分组成，伞面价格由伞半径r；绳索价格为4元/m；其它费用200元。

降落伞在降落过程中收到空气的阻力，为了确定阻力的大小，用半径3m、载重300kg的降落伞送500m高度做降落实验，测得各时刻的高度，见表2.

试确定降落伞的选购方案，即共需要多少个，每个伞的半径多大（在表1

中选择），在满足空投要求的条件下，使费用最低。

表1 不同半径的降落伞伞面价格

半径r/m 2 2.5 3 3.5 4 费用/元 65 170 350 660

1000

表

2 降落伞试验的时刻t 与高度x 的观测值 t/s 0

3 6 9 12 15 18 21 2

4 27 30 x/m 500 470 42

5 372 317 264 215 160

108

55

1

三、实验过程

模型假设

1.伞面价格1c 与伞半径r 的关系，用幂函数1b c ar =（,a b 待定）按表1数据拟合；载重m 位于球心正下方球面处，每根绳索的长度2l r =。

2.降落伞在空中只受到向下的重力和向上的空气阻力的作用，阻力与降落伞速度和伞面积的乘积成正比，阻力系数用表2数据做拟合，降落伞初速为零。 模型建立

1.目标函数

n 个降落伞的总费用记做C ，每个降落伞的费用由伞面价格1b c ar =，绳索价格

2416290.5c r r =??=和其它费用3200c =组成，于是

123()(90.5200)b C n c c c n ar r =++=++ （1）

2.伞的速度和高度

记时刻t 伞的速度为()v t ，高度为()x t ，空气阻力为2kr v ，k 是待定参数。按照牛顿第二定律，()v t 满足

2(0)0

dv m mg kr v

dt

v ?=-???=? （2）