目录

目录 .................................................................................................................................................. 1 一、Matlab基本操作与微积分计算 .............................................................................................. 4

（一）实验目的 ....................................................................................................................... 5 （二）实验学时 ....................................................................................................................... 5

1.1 求下列函数的导数. ................................................................................................... 5 1.2 求下列参数方程所确定的函数或隐函数的导数. ................................................... 5

1.3 设y?excosx,求y(4). ............................................................................................ 6 1.4 求下列多元隐函数的偏导数

?z?z,. ..................................................................... 6 ?x?yu?ln(x?a)2?(y?b)2a,b1.5 证明函数(为常数)满足拉普拉斯方程: ............. 7

1.6 计算下列积分. ........................................................................................................... 7

1?x2?x31.7求?, 并对结果求导.............................................................................. 8 241x?xt1.8 计算二重积分 ............................................................................................................ 8 1.9 计算

?L(x2?y2)dx?(x2?y2)dy,其中L为

y?1?1?x,0?x?2沿x增大

方向。 ............................................................................................................................... 9

二、 M文件、函数作图与数据拟合 ............................................................................................ 9

（一）实验目的 ....................................................................................................................... 9 （二） 实验学时 ................................................................................................................... 10

2.1画出y?secx在[0,?]之间的图象. ...................................................................... 10 2.2在同一坐标系中画出y?23x,y?x2,y?3x,y?x3,y?x的图象. ................ 10

232.3画出f(x)?(1?x)?(1?x)的图象,并根据图象特点指出函数f(x)的奇偶性......................................................................................................................................... 11 2.5已知一室模型快速静脉注射下的血药浓度数据(t=0注射300mg) ...................... 12 2.6利用教材P11页表2数据（美国人口数据，单位：百万）拟合人口指数增长模型，并检验，分析所建立的数学模型(x(t)?x0ert). ................................................ 13 2.7利用教材P29页表2数据（一盘录像带的实测数据）拟合录像带计数器的时间与计数的数学模型，并检验，分析所得数学模型(n?at?bt). ........................... 13

三、线性规划模型 ......................................................................................................................... 14

（一）实验目的 ..................................................................................................................... 14

第 1 页 共 46 页

2

（二）实验学时 ..................................................................................................................... 14

3.1 用Matlab/Lingo求解教材4.1中两个线性规划问题。 ....................................... 14 3.2 某厂生产甲乙两种口味的饮料, ............................................................................. 15 3.3 某车间有甲、乙两台机床，可用于加工三种工件 .............................................. 17 3.4 某部门在今后5年内考虑给下列项目投资， ...................................................... 18

四、非线性规划模型 ..................................................................................................................... 19

（一）实验目的 ..................................................................................................................... 19 （二）实验学时 ................................................................................................................... 19

4.1 炼油厂将A、B、C三种原料加工成甲乙丙三种汽油。 .................................. 19

五、常微分方程与级数 ................................................................................................................. 21

（一）实验目的 ..................................................................................................................... 21 （二） 实验学时 ................................................................................................................... 21

5.1 求(1)-(4)题微分方程的通解 ................................................................................... 21 5.2 求解下列初值问题 .................................................................................................. 22

xxf(x)?esinx?2cosx在点x?0的7阶taylor展开式以及在x=15.3 给出函数

处的5阶taylor展开式． ............................................................................................. 22

5.4用软件求解5.1传染病模中的相关模型。 ............................................................ 23

六、矩阵运算与层次分析模型 ..................................................................................................... 24

（一）实验目的 ..................................................................................................................... 24 （二）实验课时 ..................................................................................................................... 24

6.1 计算 .......................................................................................................................... 24

?310??102?????6.2 设A???121?,B???111?,求满足关系3A?2X?B的X．..... 25

?342??211?????6.3 求下列矩阵的秩． .................................................................................................. 25

6.4 求下列矩阵的行列式，如可逆，试用不同的方法求其逆矩阵． ...................... 26

?11?1??113?????432210?????111??125??求X． .............................................................. 27 ?=?6.5设X?6.6 解下列线性方程组． .............................................................................................. 27

6.7 （选择旅游地）某人就景色（C1），费用（C2） ............................................... 27 6.8设一容积为V(单位：m3)的大湖受到某种物质的污染， ................................... 32

七、插值与拟合 ............................................................................................................................. 33

（一）实验目的 ..................................................................................................................... 33 （二）实验课时 ..................................................................................................................... 33

7.1在化工生产中常常需要知道丙烷在各种温度T和压力P下的导热系数K． ... 33 7.2弹簧在力F的作用下伸长x，一定范围内服从胡克定律 .................................... 34

八、概率统计、回归分析模型 ..................................................................................................... 36

（一）实验目的 ..................................................................................................................... 36 （二）实验课时 ..................................................................................................................... 36

8.1 某校60名学生的一次考试成绩如下: ................................................................... 36

第 2 页 共 46 页

8.2 汽油价格 .................................................................................................................. 38 8.3根据牙膏销售量与销售差价、广告费用等数据（见下表）,建立数学模型 ...... 39

九、计算机模拟 ............................................................................................................................. 42

（一）实验目的 ..................................................................................................................... 42 （二）实验课时 ..................................................................................................................... 42

9.1.当天生产的产品必须售出，否则就会变质 ........................................................... 42 9.2.某报童以每份0.03元的价格买进报纸，以0.05元的价格出售． ..................... 43

第 3 页 共 46 页

一、Matlab基本操作与微积分计算

第 4 页 共 46 页

（一）实验目的

1.学习matlab的求导命令与求导法.

2.通过本实验加深理解积分理论中分割、近似、求和、取极限的思想方法. 3.学习并掌握用matlab求不定积分、定积分、二重积分、曲线积分的方法. 4.学习matlab命令sum、symsum与int.

（二）实验学时

1.1 求下列函数的导数.

y?(x?1)((1) （1）

1?1)22y?ln(x?x?a) y?xsinxlnxx (2) (3)

>> syms x;

>> y=( sqrt(x)+1).*(1/sqrt(x)-1); >> diff(y) ans =

(1/x^(1/2) - 1)/(2*x^(1/2)) - (x^(1/2) + 1)/(2*x^(3/2))

（2）

>> syms x;

>> y=x.*sin(x).*log(x); >> diff(y) ans =

sin(x) + log(x)*sin(x) + x*cos(x)*log(x)

（3）

>> syms a x;

>> y=log(x+sqrt(x.^2+ a.^2)); >> diff(y) ans =

(x/(a^2 + x^2)^(1/2) + 1)/(x + (a^2 + x^2)^(1/2))

1.2 求下列参数方程所确定的函数或隐函数的导数.

?x?t4?x?ln(1?t2)y22??arctg?lnx?yyxy?4ty?t?arctgtx?y??x (1) (2) (3) (4)

（1）

>> syms t;

第 5 页 共 46 页

>> x=t.^4; >> y=4*t;

>> diff(y)/diff(x) ans =

1/t^3 （2）

>> syms t;

>> x=log(1+t.^2); >> y=t-atan(t); >> diff(y)/diff(x)

ans =

-((t^2 + 1)*(1/(t^2 + 1) - 1))/(2*t)

（3）

>> syms y x;

>> f(x)=atan(y/x)-log(sqrt(x^2+y^2)); >> jacobian(f(x),[x,y])

ans =

[ - x/(x^2 + y^2) - y/(x^2*(y^2/x^2 + 1)), 1/(x*(y^2/x^2 + 1)) - y/(x^2 + y^2)]

（4）

>> syms x y; >> f(x)=x^y-y^x;

>> jacobian(f(x),[x,y])

ans =

[ x^(y - 1)*y - y^x*log(y), x^y*log(x) - x*y^(x - 1)]

1.3 设y?excosx,求y(4).

>> syms x;

>> y=exp(x)*cos(x); >> diff(y,x,4)

ans =

-4*exp(x)*cos(x)

1.4 求下列多元隐函数的偏导数

?z?z,. ?x?y222zcosx?cosy?cosz?1e(1) (2)?xyz

第 6 页 共 46 页

（1）

>> syms x y z;

>> f(x)=cos(x)^2+cos(y)^2+cos(z)^2-1; >> a=[diff(f(x),x),diff(f(x),y),diff(f(x),z)]; >> dz_dx=a(1)/a(3)

d z_dx =

(cos(x)*sin(x))/(cos(z)*sin(z))

>> dz_dy=a(2)/a(3)

dz_dy =

(cos(y)*sin(y))/(cos(z)*sin(z)) a =

[ -2*cos(x)*sin(x), -2*cos(y)*sin(y), -2*cos(z)*sin(z)] （2）

>> syms x y z;

>> f(x)=exp(z)-x*y*z;

>> a=[diff(f(x),x),diff(f(x),y),diff(f(x),z)]; >> dz_dx=a(1)/a(3)

dz_dx =

-(y*z)/(exp(z) - x*y)

>> dz_dy=a(2)/a(3)

dz_dy =

-(x*z)/(exp(z) - x*y)

u?ln(x?a)2?(y?b)2a,b1.5 证明函数(为常数)满足拉普拉斯方程:

?2u?2u?2?02?x?y （提示:对结果用simplify化简）

syms x y a b

du_dx=diff(log(sqrt((x-a)^2+(y-b)^2)),x,2); du_dy=diff(log(sqrt((x-a)^2+(y-b)^2)),y,2); simplify(du_dx+du_dy)

ans = 0

1.6 计算下列积分.

(1)

?sin2xdx1?sin2x1 （2）

?3dx(2x2?1)1?x2?0 (3)

??34xdx sin2x第 7 页 共 46 页

（1）

>> syms x;

>> syms a real;

>> y=sin(2*x)/(sqrt(1+sin(x).^2)); >> int(y,x) ans =

2*(sin(x)^2 + 1)^(1/2) （2）

>> syms x;

>> y=1/(2*(x^2)+1)*(sqrt(1+x^2)); >> int(y,0,(1/sqrt(3))) ans =

pi/4 + asinh(3^(1/2)/3)/2 - atan(((2*2^(1/2))/5 + 3^(1/2)/10))/2

（3）

>> syms x;

>> y=x/(sin(x))^2; >> int(y,pi/4,pi/3) ans =

pi/4 + log((2^(1/2)*3^(1/2))/2) - (pi*3^(1/2))/9

?t1?x2?x31.7求1x2?x4, 并对结果求导

>> syms x; >> syms t real;

>> y=(1+x^2+x^3)/(x^2+x^4); >> a=int(y,1,t) a =

log(t^2 + 1)/2 - log(2)/2 - 1/t + 1

>> diff(a)

ans =

t/(t^2 + 1) + 1/t^2

1.8 计算二重积分

??(x?y)dxdy??(x2?y2)dxdy(1)

x2?y2?1 (2)

x2?y2?x

第 8 页 共 46 页

3^(1/2)/5)/(2^(1/2)/5 +

(3)

??xydxdyD2y?x?1y?2x?6所围成的闭区域。 ,其中D是由直线与抛物线

xsindxdy??3yy?x,y?2和曲线x?yD(4) ,其中D是由直线所围成闭区域。

（1） >> syms x >> syms y

>> int(int((x+y),y,-sqrt(1-x^2),sqrt(1-x^2)),x,-1,1) ans = 0 （2）

>> syms x; >> syms y;

>> int(int((x^2+y^2),y,-sqrt(x-x^2),sqrt(x-x^2)),x,0,1) ans =

(3*pi)/32 （3） syms x y

int(int(x*y,y,-sqrt(2*x+6),sqrt(2*x+6)),x,-3,-1) +int(int(x*y,y,x-1,sqrt(2*x+6)),x,-1,36 (4) syms x y

int(int(sin(x/y),y,x^(1/3),x),x,1,2)+int(int(sin(x/y),y,x^(1/3),2),x,2,8) (3*cos(1))/2 + sin(1)/2 - sin(4)/2

1.9 计算

?L(x2?y2)dx?(x2?y2)dy,其中L为

y?1?1?x,0?x?2沿x增大方向。

syms x

int(2*x^2,0,1)+int(2*x^2-8*x+8,1,2) 4/3

二、 M文件、函数作图与数据拟合

（一）实验目的

1学习matlab自定义函数.

2加深理解罗必塔法则、极值、最值、单调性.

3学习matlab函数绘图命令(二维、三维图形，散点图). 学会使用matlab帮助系统。

第 9 页 共 46 页

4掌握matlab基本数据拟合命令，用数据拟合解决一些实际问题.

（二）实验学时

2.1画出y?secx在[0,?]之间的图象.

>> x=0:pi; >> y=sec(x);

>> ezplot('sec(x)')

2.2在同一坐标系中画出y?x,y?x2,y?3x,y?x3,y?x的图象.

x=1:0.01:100;y1=x.^(1/2);hold on; >> y3=x.^(1/3); >> y4=x.^3; >> y5=x; >> y2=x.^2;

>> plot(x,y1,'r',x,y2,'b',x,y3,'g',x,y4,'c',x,y5,'y')

第 10 页 共 46 页

2323

2.3画出f(x)?(1?x)?(1?x)的图象,并根据图象特点指出函数f(x)的奇偶性

x=-10:1:10;

>> y=(1-x).^(2/3)+(1+x).^(2/3); >>plot(x,y);

2.4已知热敏电阻数据：

第 11 页 共 46 页

温度t(0C) 电阻R(?) 20.5 765 32.7 826 51.0 873 73.0 942 95.7 1032 求6000C时的电阻R。(设R=at+b,a,b为待定系数)

解：

>> syms t;

x=[20.5,32.7,51.0,73.0,95.7]; y=[765,826,873,942,1032]; A=polyfit(x,y,1) a=[t;1]; R=A*a b=[600;1]; R_600=A*b

A =

3.3987 702.0968

R =

(1913320152740319*t)/562949953421312 + 3087854411876483/4398046511104 R_600 =

2.7413e+03

2.5已知一室模型快速静脉注射下的血药浓度数据(t=0注射300mg) t (h) c(t) (?g/ml) 0.25 19.21 0.5 18.15 1 15.36 1.5 14.10 2 12.89 3 9.32 4 7.45 6 5.24 8 3.01 求血药浓度随时间的变化规律c(t). (提示：t 与y=ln(c(t))呈线性关系)

解：

syms t;

format long;

x=[0.25,0.5,1,1.5,2,3,4,6,8];

y=[19.21,18.15,15.36,14.10,12.89,9.32,7.45,5.24,3.01]; f=log(y);

A=polyfit(x,f,1) a=[t;1]; R=A*a; c_t=exp(R) >>A =

-0.234718197179491 2.994276199298933

c_t =

exp(3371255293851755/1125899906842624-(4228307141418513*t)/18014398509481984)

第 12 页 共 46 页

2.6利用教材P11页表2数据（美国人口数据，单位：百万）拟合人口指数增长模型，并检验，分析所建立的数学模型(x(t)?x0ert). 年份 实际人口 年份 1790 3.9 1900 1800 5.3 1910 1810 1820 1830 01840 7.2 9.6 12.9 17.1 1850 1860 1870 1880 1890 23.2 31.4 38.6 50.2 62.9 1920 1930 1940 1950 1960 1970 198 1990 2000 实际人76.0 92. 106.123.131.150.7 179.204.226.251.281.5 2 7 3 0 5 4 4 口 解： syms t;

format long; x=1790:10:2000;

y=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6,50.2,62.9,76.0,92.0,106.5,123.2,131.7,150.7,1

79.3,204.0,226.5,251.4,281.4];

f=log(y);

B=polyfit(x,f,4) a=[t^4;t^3;t^2;t;1]; R=B*a; x_t=B*a

z=polyval(B,x); plot(x,f,'k+',x,z,'r') err=norm(z-f,2) >>B =

1.0e+004 *

0.000000000000292 -0.000000002205327 0.000006241457716

-0.007836484321485 3.682281541453110

x_t =

(3527067531156703*t^4)/1208925819614629174706176 - (1627243798949389*t^3)/73786976294838206464 + (8994888525520685*t^2)/144115188075855872 - (689304450588591*t)/8796093022208 + 632611151808223/17179869184

err =

0.103749826591387

2.7利用教材P29页表2数据（一盘录像带的实测数据）拟合录像带计数器的时间与计数的

2n?at?bt). 数学模型，并检验，分析所得数学模型(

t(分) n t(分) n 0 0000 100 4004 10 0617 110 4280 20 1141 120 4545 30 1601 130 4803 40 2019 140 5051 50 2403 150 5291 60 2760 160 5525 70 3096 170 5752 80 3413 184 6061 90 3715 第 13 页 共 46 页

三、线性规划模型

（一）实验目的

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

2、掌握用数学软件(Matlab/Lingo)求解线性规划问题;

（二）实验学时

3.1 用Matlab/Lingo求解教材4.1中两个线性规划问题。 （1）

maxz?72x1?64x2s.t.x1?x2?5012x1?8x2?4803x1?100x1,x2?0（2）

maxz?24x1?16x2?44x3?32x4?3x5?3x6s.t.x1?x5x2?x6??50344(x1?x5)?2(x2?x6)?2x5?2x6?480x1?x5?100x3?0.8x5x4?0.75x6

x1,x2,x3,x4,x5,x6?0解： (1)

format long; c=[-72 -64];

A=[1 1;12 8;3 0]; b=[50;480;100]; Aeq=[]; beq=[]; vlb=[0;0]; vub=[];

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

Optimization terminated.

x =

19.999999996748322

第 14 页 共 46 页

30.000000003205681

fval =

-3.359999999971043e+003 (2)

c=[-24 -16 -44 -32 3 3]; A=[4 2 0 0 6 4;1 0 0 0 1 0]; b=[480;100];

Aeq=[0 0 1 0 -0.8 0;0 0 0 1 0 -0.75]; beq=[0;0];

vlb=[0;0;0;0;0;0]; vub=[ ];

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

Optimization terminated.

x =

1.0e+002 *

0.000000000000067 2.399999999999335 0.000000000000014 0.000000000000157 0.000000000000017 0.000000000000210

fval =

-3.839999999999591e+003

3.2 某厂生产甲乙两种口味的饮料,每百箱甲饮料需用原料6千克,工人10名,可获利10万元;每百箱乙饮料需用原料5千克,工人20名,可获利9万元.今工厂共有原料60千克,工人150名,又由于其他条件所限甲饮料产量不超过8百箱.问如何安排生产计划,即两种饮料各生产多少使获利最大.进一步讨论:

1)若投资0.8万元可增加原料1千克,问应否作这项投资. 2)若每百箱甲饮料获利可增加1万元,问应否改变生产计划. max z?10x?9y

6x?5y?6010x?20y?150x?8 s.t. x,y?0

解：

format long; c=[-10 -9];

A=[6 5;10 20;1 0]; b=[60;150;8]; Aeq=[]; beq=[]; vlb=[0;0];

第 15 页 共 46 页

vub=[];

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

Optimization terminated.

x =

6.428571428446904 4.285714285701694

fval =

-1.028571428557843e+002 甲：6 乙：4

总获利：96 （1）

max z?10x?9y

6x?5y?6110x?20y?150x?8s.t. x,y?0

c=[-10 -9];

A=[6 5;10 20;1 0]; b=[61;150;8]; Aeq=[]; beq=[]; vlb=[0;0]; vub=[];

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

Optimization terminated. x =

6.714285714165383 4.142857142846659

fval =

-1.044285714272738e+002 甲：6 乙：4

总获利：96-0.8=95.2 （2)

max z?11x?9y

6x?5y?6010x?20y?150x?8s.t. x,y?0

第 16 页 共 46 页

c=[-11 -9];

A=[6 5;10 20;1 0]; b=[60;150;8]; Aeq=[]; beq=[]; vlb=[0;0]; vub=[];

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

Optimization terminated. x =

7.999999999630490 2.400000000370391

fval =

-1.095999999992689e+002 甲：8 乙：2

总获利：106

所以改变生产计划。

3.3 某车间有甲、乙两台机床，可用于加工三种工件。假定这两台车床的可用台时数分别为800和900，三种工件的数量分别为400、600和500，且已知用三种不同车床加工单位数量不同工件所需的台时数和加工费用如下表。问怎样分配车床的加工任务，才能既满足加工工件的要求，又使加工费用最低？

单位工件所需加工台时数 车床类型 工件1 甲 乙 min

0.4 0.5 工件2 1.1 1.2 工件3 1.0 1.3 工件1 13 11 工件2 9 12 工件3 10 8 800 900 单位工件的加工费用 可用台时数 z?13x1?9x2?10x3?11x4?12x5?8x6

0.4x1?1.1x2?x3?8000.5x4?1.2x5?1.3x6?900x1?x4?400x2?x5?600x3?x6?500s.t

x1,x2,x3,x4,x5,x6?0

解：

c=[13 9 10 11 12 8];

A=[0.4 1.1 1 0 0 0;0 0 0 0.5 1.2 1.3;-1 0 0 -1 0 0;0 -1 0 0 -1 0;0 0 -1 0 0 -1]; b=[800;900;-400;-600;-500];

第 17 页 共 46 页

Aeq=[]; beq=[];

vlb=[0;0;0;0;0;0]; vub=[];

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

Optimization terminated. x =

0.0000 600.0000 0.0000 400.0000 0.0000 500.0000

fval =

1.3800e+004

分配车床的加工任务如下：

单位工件所需加工台时数 车床类型 工件1 甲 乙 0.4 0.5 工件2 1.1 1.2 工件3 1.0 1.3 工件1 13 11 工件2 9 12 工件3 10 8 800 900 单位工件的加工费用 可用台时数 才能既满足加工工件的要求，又使加工费用最低且最低费用是13800。 3.4 某部门在今后5年内考虑给下列项目投资，已知：

项目A：从第一年到第四年初投资，并于次年底回收本例115%；

项目B：从第三年初投资，到第五年底回收本利125%，投资额不超过4万元。 项目C：从第二年初投资，到第五年底回收本利140%，投资额不超过3万元。 项目D：五年内每年初购公债，于本年底归还本利，年利为6%。

该本部门现有资金10万元，问应如何确定每年各项目投资额，使到第五年底，部门拥有资金的总额最大。

解：

c=[-1.3225 -1.25 -1.4 -1.33822558]; A=[1 1 1 1;0 1 0 0;0 0 1 0]; b=[10;4;3]; Aeq=[]; beq=[];

vlb=[0;0;0;0]; vub=[];

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

Optimization terminated. x =

第 18 页 共 46 页

0.0000 0.0000 3.0000 7.0000

fval =

-13.5676

四、非线性规划模型

（一）实验目的

1、直观了解非线性规划的基本内容。 2、掌握用数学软件求解优化问题。

（二）实验学时

4.1 炼油厂将A、B、C三种原料加工成甲乙丙三种汽油。一桶原油加工成汽油的费用为4元，每天至多能加工汽油14，000桶。原油的买入价、买入量、辛烷值、硫含量，及汽油的卖出价、需求量、辛烷值、硫含量由表4-15给出。问如何安排生产计划，在满足需求的条件下使利润最大？

一般来说，作广告可以增加销售，估计一天向一种汽油投入一元广告费，可使该汽油日销售量增加10桶，且每天最多投入广告费800元，问：如何安排生产和广告计划使利润最大？

一． 问题重述

已知某炼油厂原材料和产品的价格，辛烷值，含量及数量的约束条件，要求在一定

的加工费用下求出如何购进原材料，如何分配原材料才能使这个工厂的利润最大。

(1)由上表给出，如何安排生产计划，使利润最大？ (2)投入广告后，如何安排生产和广告计划使利润最大？ 二．符号说明

f(x)：未投入广告前的最大利润； F(x)：投入广告后的最大利润； M：卖出汽油获得的总收入； N：买入原油需要的总费用； P：原油加工成汽油的总费用； Q：投入的总的广告费；

A1、A2、A3：分别是由A类原油加工成甲、乙、丙三种汽油的桶数； B1、B2、B3：

分别是由B类原油加工成甲、乙、丙三种汽油的桶数； C1、C2、C3：分别是由C类原油加工成甲、乙、丙三种汽油的桶数； X1、X2、X3：分别是投入到甲、乙、丙三种汽油上的广告费用。

三． 建立模型及模型的求解 （1）未投入广告前

卖出汽油获得的总收入：M=70（A1+B1+C1）+60（A2+B2+C2）+50（A3+B3+C3）。 买入原油需要的总费用：N=45（A1+A2+A3）+35（B1+B2+B3）+25（C1+C2+C3）。

原油加工成汽油的总费用：P=4（A1+A2+A3+ B1+B2+B3+ C1+C2+C3）。

约束条件：

1. 由题知每天至多能加工汽油14，000桶，故：

第 19 页 共 46 页

A1+A2+A3+B1+B2+B3+ C1+C2+C3<=14000。 2. 规定每天每种原油最多买入5000桶，故： A1+A2+A3<=5000； B1+B2+B3<=5000； C1+C2+C3<=5000。

3. 关于需求量分为两种情况：(1)产量大于需求量；(2)产量等于需求量。 对于第一

种情况，约束条件为：

A1+B1+C1>3000； A2+B2+C2>2000； A3+B3+C3>1000。

对于第二种情况，约束条件为： A1+B1+C1=3000； A2+B2+C2=2000； A3+B3+C3=1000。

4. 生产出三种汽油对辛烷值有如下要求： 12A1+6B1+8C1>=10(A1+B1+C1)； 12A2+6B2+8C2>=8(A2+B2+C2)； 12A3+6B3+8C3>=10(A3+B3+C3)。

5. 生产出三种汽油对硫含量有如下要求： 0.5A1+2B1+3C1<=A1+B1+C1； 0.5A2+2B2+3C2<= A2+B2+C2； 0.5A3+2B3+3C3<= A3+B3+C3。

目标函数：未投入广告前的最大利润=卖出汽油获得的总收入-买入原油需要的总费

用-原油加工成汽油的总费用;

即f(x)=M-N-P=70（A1+B1+C1）+60（A2+B2+C2）+50（A3+B3+C3）-45（A1+A2+A3）

+35（B1+B2+B3）+25（C1+C2+C3）-4（A1+A2+A3+ B1+B2+B3+ C1+C2+C3）= 21A1+11A2+A3+ 31B1+21B2+11B3+ 41C1+31C2+21C3

☉用lingo对模型求解,结果如下：

(1)投入广告前产量大于需求量时的最大利润为：142500。

此时每天的生产计划为：由A类原油加工成甲、乙、丙三种汽油的量为3000、1333、

666；由B类原油加工成甲、乙、丙三种汽油的量为：1500、667、333；不由C类原油加工成甲、乙、丙三种汽油。

(2)投入广告前产量等于需求量时的最大利润为：110000。

此时每天的生产计划为：由A类原油加工成甲、乙、丙三种汽油的量为2400、1600、

800；不由B类原油加工成甲、乙、丙三种汽油；由C类原油加工成甲、乙、丙三种汽油的量为600、400、200。

（2）投入广告后

投入的总的广告费：Q=X1+X2+X3。 此时约束条件3发生改变：

产量大于需求量时的约束条件为： A1+B1+C1>3000+10X1； A2+B2+C2>2000+10X2； A3+B3+C3>1000+10X3。

产量等于需求量时的约束条件为： A1+B1+C1=3000+10X1； A2+B2+C2=2000+10X2； A3+B3+C3=1000+10X3。

目标函数：投入广告后的最大利润=卖出汽油获得的总收入-买入原油需要的总费用

-原油加工成汽油的总费用-投入的总的广告费;

即F(x)=M-N-P-Q=70（A1+B1+C1）+60（A2+B2+C2）+50（A3+B3+C3）-45

（A1+A2+A3）+35（B1+B2+B3）+25（C1+C2+C3）-4（A1+A2+A3+ B1+B2+B3+ C1+C2+C3）

第 20 页 共 46 页

-X1+X2+X3= 21A1+11A2+A3+31B1+21B2+11B3+ 41C1+31C2+21C3 -X1+X2+X3

☉用lingo对模型求解，结果如下：

(1)投入广告后产量大于需求量时的最大利润为：142500。

此时每天的生产计划为：由A类原油加工成甲、乙、丙三种汽油的量为3000、1333、

667；由B类原油加工成甲、乙、丙三种汽油的量为1500、667、333；不由C类原油加工成甲、乙、丙三种汽油。

不投入广告到甲、乙、丙三种汽油上。

(2)投入广告后产量等于需求量时的最大利润为：142350

此时每天的生产计划为：由A类原油加工成甲、乙、丙三种汽油的量为3000、1333、

667；由B类原油加工成甲、乙、丙三种汽油的量为1500、667、333；不由C类原油加工成甲、乙、丙三种汽油。

分别投入到甲、乙、丙三种汽油上的广告费用为150、0、0。 四． 结果分析

由上可见，投入广告前最大利润为：110000。投入广告后最大利润为：142350。故

投入广告后能为工厂增加更多利润。而且为了达到最大利润，投入广告前不由B类原油加工成甲、乙、丙三种汽油，不由C类原油加工成甲、乙、丙三种汽油，即可以不用某种原油。

五、常微分方程与级数

（一）实验目的

１．学习用matlab求解微分方程命令dsolve. ２．学习matlab泰勒级数展开命令. ３．巩固幂级数的收敛半径、和等概念.

（二）实验学时

5.1 求(1)-(4)题微分方程的通解

22(1) 2xyy??y?1 (2) (xcosy?sin2y)y??1 x(3) y???3y??y?ecos2x (4) y???4y?x?1?sinx

解： (1)

y=dsolve('Dy*2*(x^2)*y=y^2+1','x') y =

i -i (exp(C1 - 1/x) - 1)^(1/2) -(exp(C1 - 1/x) - 1)^(1/2) (2)

dsolve('(x*cos(y)+sin(2*y))*Dy=1','x')

第 21 页 共 46 页

ans =

-asin(lambertw(-1/2*C1*exp(-1/2*x-1))+1/2*x+1) (3)

dsolve('D2y+3*Dy-y=exp(x)*cos(2*x)','x')

ans =

C2*exp(x*(13^(1/2)/2 - 3/2)) + C3/exp(x*(13^(1/2)/2 + 3/2)) + (13^(1/2)*exp((5*x)/2 -

(13^(1/2)*x)/2)*exp(x*(13^(1/2)/2 - 3/2))*(2*sin(2*x) - cos(2*x)*(13^(1/2)/2 - 5/2)))/(13*((13^(1/2)/2 - 5/2)^2 + 4)) - (13^(1/2)*exp((5*x)/2 + (13^(1/2)*x)/2)*(2*sin(2*x) + cos(2*x)*(13^(1/2)/2 + 5/2)))/(13*exp(x*(13^(1/2)/2 + 3/2))*((13^(1/2)/2 + 5/2)^2 + 4))

(4)

dsolve('D2y+4*y=x+1+sin(x)','x')

ans =

cos(2*x)*(cos(2*x)/4 - sin(2*x)/8 + sin(3*x)/12 - sin(x)/4 + (x*cos(2*x))/4 - 1/4) +

sin(2*x)*(cos(2*x)/8 - cos(3*x)/12 + sin(2*x)/4 + cos(x)/4 + (x*sin(2*x))/4 + 1/8) + C5*cos(2*x) + C6*sin(2*x) 5.2 求解下列初值问题

?d2xdx?2n?a2x?0?2?2dt222dy?dt?x?2xy?y?(y?2xy?x)?0?dx??x?x,dx?Vt?00t?00?yx?1?1?dt??(1) (2)

解：

(1)

dsolve('x^2+2*x*y-y^2+(y^2+2*x*y-x^2)*Dy=0','y(1)=1','x') ans =

x*((4/x + 1/x^2 - 4)^(1/2)/2 + 1/(2*x)) (2)

dsolve('D2y+2*n*Dy+(a^2)*y=0','y(0)=x_0','Dy(0)=V','t') ans =

(V + n*x_0 + x_0*(n^2 - a^2)^(1/2))/(2*exp(t*(n - (n^2 - a^2)^(1/2)))*(n^2 -

a^2)^(1/2)) - (V + n*x_0 - x_0*(n^2 - a^2)^(1/2))/(2*exp(t*(n + (n^2 - a^2)^(1/2)))*(n^2 - a^2)^(1/2))

xxf(x)?esinx?2cosx在点x?0的7阶taylor展开式以及在x=1处的55.3 给出函数

阶taylor展开式．

解： syms x;

f=exp(x)*(sin(x))+(cos(x))*(2^x); taylor(f,0,8) taylor(f,1,6) ans =

(log(2)^3/144 - log(2)/720 - log(2)^5/240 + log(2)^7/5040 - 1/630)*x^7 + (log(2)^2/48 -

第 22 页 共 46 页

log(2)^4/48 + log(2)^6/720 - 1/80)*x^6 + (log(2)/24 - log(2)^3/12 + log(2)^5/120 - 1/30)*x^5 + (log(2)^4/24 - log(2)^2/4 + 1/24)*x^4 + (log(2)^3/6 - log(2)/2 + 1/3)*x^3 + (log(2)^2/2 + 1/2)*x^2 + (log(2) + 1)*x + 1

ans =

64*cos(6) + exp(6)*sin(6) 5.4用软件求解5.1传染病模中的相关模型。 一、SI模型 假设

1) 总人数N不变，t 时刻病人和健康人的比例分别为i(t),s(t); 2) 每个病人每天有效接触人数为?，且使接触的健康人致病 则可以建立如下模型

?di???i(1?i) ?dt??i(0)?i0解：

0 i=dsolve('Dy=a*y*(1-y)','y(0)=b','t')%

i =

-1/(exp(log((b - 1)/b) - a*t) - 1)

二、SIS模型

若传染病无免疫性——病人治愈成为健康人，健康人可再次被感染。 增加假设3）病人每天治愈的比例为?。

a??,b?i则可以建立如下模型

?di???i(1?i)??i ?dt??i(0)?i0???/? ~ 一个感染期内每个病人的有效接触人数，称为接触数，即 1?di???i[i?(1?)]?? ?dt??i(0)?i0 a??,u??,b?i0

解：

i=dsolve('Dy=a*y*(1-y)-u*y','y(0)=b','t')% i =

-((tan((t - (2*atan((- a*i + u*i + 2*a*b*i)/(a - u))*i)/(a - u))*((a*i)/2 - (u*i)/2)) + (a*i -

u*i)/(a - u))*(a - u)*i)/(2*a)

第 23 页 共 46 页

六、矩阵运算与层次分析模型

（一）实验目的

1理解矩阵及数组概念．

2掌握matlab对矩阵及数组的操作命令（加法、数乘、乘法、转置、逆）． 3求AX=B的通解．

4掌握matlab求矩阵的秩命令． 5掌握matlab求方阵的行列式命令．

6理解逆矩阵概念，掌握matlab求逆矩阵命令． 7会用matlab求解线性方程组．

（二）实验课时

6.1 计算

?123??214???1??02?10?12?????1?12?2?023??＋??； （1）??105?????10?020?????312?1??101??15???02????0310030?； ????（2）?

5?123??214???1??02?10?12?????1?12?2?023??＋??； （1）?

A=[1,2,3;0,2,-1;1,-1,2]; >> B=[2,1,4;0,-1,2;0,2,3]; >> A+1/2*(B.^5)

ans =

17.0000 2.5000 515.0000 0 1.5000 15.0000 1.0000 15.0000 123.5000

5第 24 页 共 46 页

?1??0?312?1??1???0310???0（2）

02035??0?1??0???10???15???02???；

A=[3,1,2,-1;0,3,1,0];

>> B=[1,0,5;0,2,0;1,0,1;0,3,0]; >> C=[-1,0;1,5;0,2]; >> A*B*C

ans =

-6 29 5 32

?310??102?????6.2 设A???121?,B???111?,求满足关系3A?2X?B的X．

?342??211?????

A=[3,1,0;-1,2,1;3,4,2]; >> B=[1,0,2;-1,1,1; 2,1,1]; >> (3*A+B)/2

ans =

5.0000 1.5000 1.0000 -2.0000 3.5000 2.0000 5.5000 6.5000 3.5000

6.3 求下列矩阵的秩．

?120??? (1) ?011? (2)

??123????25??75?75??25319494321743??53132?

54134??2048?1）a=[1,2,0;0,1,1;-1,2,3]; rank(a)

ans =

3

2）b=[25,31,17,43;75,94,53,132;75,94,54,134;25,32,20,48]; rank(b)

第 25 页 共 46 页

ans = 3

6.4 求下列矩阵的行列式，如可逆，试用不同的方法求其逆矩阵．

?22?1???（1）?1?24?

?582????1??2?1?1 （2）?4??312?11?1??0?2?6? 231）a=[2,2,-1;1,-2,4;5,8,2]; >> det(a)

ans =

-54

>> inv(a)

ans =

0.6667 0.2222 -0.1111 -0.3333 -0.1667 0.1667 -0.3333 0.1111 0.1111 2）

a=[1,2,3,4;2,3,1,2;1,1,1,-1;1,0,-2,-6]; >> det(a)

ans =

-1.0000

>> inv(a)

ans =

22.0000 -6.0000 -26.0000 17.0000 -17.0000 5.0000 20.0000 -13.0000 -1.0000 -0.0000 2.0000 -1.0000 4.0000 -1.0000 -5.0000 3.0000

第 26 页 共 46 页

?11?1??113?????

432210????

?111??125?

?求X． ?=?6.5设X?a=[1,1,-1;2,1,0;1,1,1];

>> b=[1,1,3;4,3,2;1,2,5]; >> b*inv(a)

ans =

-1 0 2 0 1 2 -1 -1 4

6.6 解下列线性方程组．

?2x1?2x2?x3?x4?4?4x?3x?x?2x?6?1234??8x1?3x2?3x3?4x4?12?3x?3x2?2x2?2x4?6 (1) ?1

?x1?2x2?x3?x4?x5?0?2x?x?x?2x?3x?0?12345 (2) ?

?3x1?2x2?x3?x4?2x5?0??2x1?5x2?x3?2x4?2x5?0（1）

a=[2,2,-1,1;4,3,-1,2;8,3,-3,4;3,3,-2,-2]; >> b=[4,6,12,6]; >> x=b/a x =

-15.4286 17.2857 -3.8571 -1.1429

（2）a=[1,-2,1,-1,1,0;2,1,-1,2,-3,0;3,-2,-1,1,-2,0;2,-5,1,-2,2,0]; >> rref(a)

ans =

1.0000 0 0 0.5000 -0.8750 0 0 1.0000 0 0.5000 -0.6250 0 0 0 1.0000 -0.5000 0.6250 0 0 0 0 0 0 0

6.7 （选择旅游地）某人就景色（C1），费用（C2），居住（C3），饮食（C4），旅途（C5）五个准则，从三个旅游地桂林（P1），黄山（P2），苏州（P3）中选择最佳旅游地。其用层次分析法（AHP）进行选择，过程如下： (1)建立的成对比较矩阵

第 27 页 共 46 页

433??112??21755??11?11A??14723?11?211?35??11311??35?

(2) 用和法近似计算A的最大特征根?

a.将矩阵A按列向量归一化，得矩阵B； b.将矩阵B按行求和，得向量W；

c.将向量W归一化，得向量w，作为最大特征根的近似特征向量； d.根据Aw??w,求出?作为最大特征根的近似值； (3) 检验矩阵A的一致性

??na.计算A的一致性指标 CI=n?1

b.检验一致性比率 CR=CI /RI是否小于0.1。（其中n=5时，RI＝1.12） (4) 若通过一致性检验，将w作为准则Ci对目标层O的权向量。 (5) 同样，构造第3层(方案层)对第2层(准则层)每一准则的成对比较阵，

?1?B1??1/?1/?2?513??11/31/8??1?????21B2?311/3B?11323?????????51/2?8131?1/31/31??，??，，

?134??111/4?????B4??1/311?B5??111/4??1/411??441???，??

(6) 由

Bk计算出方案对第k个准则的权向量wk,最大特征根?k和一致性指标CIk，并

作一次性检验。（其中n=3时，RI＝0.58） (7) 计算方案对目标的组合权向量,

w3?W3w，其中W3?[w1,w2,?,w5]

然后通过组合权向量进行选择。

（建议用M文件建立用和法求权向量、最大特征根及一致性检验的函数，方便调用） （1）（2）

解：

A=[1 1/2 4 3 3;2 1 7 5 5;1/4 1/7 1 1/2 1/3;1/3 1/3 2 1 1;1/3 1/3 3 1 1]; B=A;

[m,n]=size(A); for i=1:n

B(:,i)=A(:,i)/sum(A(:,i)); end

第 28 页 共 46 页

B

W=sum(B'); W=W'

w=W/sum(W) d=A*w; d./w B =

0.2553 0.2165 0.2353 0.2857 0.2903 0.5106 0.4330 0.4118 0.4762 0.4839 0.0638 0.0619 0.0588 0.0476 0.0323 0.0851 0.1443 0.1176 0.0952 0.0968 0.0851 0.1443 0.1765 0.0952 0.0968 W =

1.2831 2.3155 0.2644 0.5391 0.5979 w =

0.2566 0.4631 0.0529 0.1078 0.1196

ans =

5.3848 5.3629 5.2378 5.3150 5.2343

（3）

A=[1 1/2 4 3 3;2 1 7 5 5;1/4 1/7 1 1/2 1/3;1/3 1/3 2 1 1;1/3 1/3 3 1 1]; [m,n]=size(A); RI=1.12; for i=1:n

B(:,i)=A(:,i)/sum(A(:,i)); end

W=sum(B'); W=W';

w=W/sum(W); d=A*w; R=d./w;

r=(R(1)+R(2)+R(3)+R(4)+R(5))/5; disp(['最大特征根 r = ',num2str(r)]) CI=(r-5)/4;

第 29 页 共 46 页

disp(['一致性指标 CI = ',num2str(CI)]) CR=CI/RI;

disp(['一致性比率 CR = ',num2str(CR)])

结果：

最大特征根 r = 5.307 %r?? 一致性指标 CI = 0.076743 一致性比率 CR = 0.068521 可知：CR<0.1 （4）w =

0.2566 0.4631 0.0529 0.1078 (5) w3 =

0.3014 0.2438 0.4548 (6)

主函数如下： function main() RI=0.58; for i=1:1:5 if i==1

A=[1 2 5;1/2 1 2;1/5 1/2 1]; [w]=SUM(A,RI) end if i==2

A=[1 1/3 1/8;3 1 1/3;8 3 1]; [w]=SUM(A,RI) end if i==3

A=[1 1 3;1 1 3;1/3 1/3 1]; [w]=SUM(A,RI) end if i==4

A=[1 3 4;1/3 1 1;1/4 1 1]; [w]=SUM(A,RI) end if i==5

A=[1 1 1/4;1 1 1/4;4 4 1]; [w]=SUM(A,RI) end end

子函数如下：

function [w]=SUM(A,RI) [m,n]=size(A); for i=1:n

B(:,i)=A(:,i)/sum(A(:,i)); end

W=sum(B'); W=W';

w=W/sum(W); d=A*w; R=d./w;

r=(R(1)+R(2)+R(3))/3;

第 30 页 共 46 页0.1196

disp(['最大特征根 r =',num2str(r)]) CI=(r-3)/2; CR=CI/RI;

disp(['一致性比率 CR =',num2str(CR)])

最大特征根 r =3.0055

一致性比率 CR =0.0047747 w =

0.5949 0.2766 0.1285

最大特征根 r =3.0015

一致性比率 CR =0.0013294 w =

0.0820 0.2364 0.6816

最大特征根 r =3 一致性比率 CR =0 w =

0.4286 0.4286 0.1429

最大特征根 r =3.0092

一致性比率 CR =0.0079433 w =

0.6327 0.1924 0.1749

最大特征根 r =3 一致性比率 CR =0 w =

0.1667 0.1667 0.6667 (7)

w_1=[0.5949;0.2766;0.1285]; w_2=[0.0820;0.2364;0.6816]; w_3=[0.4286;0.4286;0.1429]; w_4=[0.6327;0.1924;0.1749]; w_5=[0.1667;0.1667;0.6667]; W_3=[w_1,w_2,w_3,w_4,w_5]

第 31 页 共 46 页

w =[0.2566;0.4631;0.0529;0.1078;0.1196] w3=W_3*w w3 =

0.3014 0.2438 0.4548

结果表明方案P3在旅游地选择中占的权重近于1/2，远大于P1，P2，应作为第一选

择地点。

6.8设一容积为V(单位：m3)的大湖受到某种物质的污染，污染物均匀地分布在湖中．若从某时刻起污染源被切断，设湖水更新的速率是r(单位：m3/d)．试建立求污染物浓度下降至原来的5%需要多少时间的数学模型．美国密西根湖的容积为4871×109m3，湖水的流量为3.663959132×1010m3/d，求污染中止后，污染物浓度下降到原来的5%所需要的时间． 一、模型假设 1、，湖水流量为常量，湖水体积为常量； 2、流入流出湖水水污染浓度为常量， 二、符号说明

W（t）: t 时刻水污染浓度 t: 时间，以天作单位 m: 外进湖中水污染浓度 r: 湖水的更新速率 V: 湖水的体积 三、问题分析

问题（一）要求经过500天湖水的浓度，由于流入和流出的湖水浓度不同，我们在考虑此问题时，运用微积分方程和质量守恒定律得出水污染浓度与已知量之间的关系；问题（二）污染源被切断的情况，即湖水的污染浓度不再改变，即m=0，由于问题（二）给出污染物浓度下降到原来的5%，从而可以求得所需的时间。

四、模型的建立与求解

设t时刻湖区的污染物浓度为W（t），考虑时间区间[t，t +⊿t]并利用质量守恒定律： [t，t +⊿t]内湖中污染浓度的变化量=流入湖水的污染量—流出湖水的污染量。 用数学表达式表示为：

V[w(t+⊿t )—w（t）]=rm*⊿t—

于是得，令⊿t→0

dw/dt=a-bx t>0,w(0)=wo 其中a=rm,b=r

求得w(t)=m+(wo-m)e^(-r*t)……（1） 问题（一）中v=5.941e12(m3)

r=4.45365e10, wo=10%, m=5%,t=500代入式中，得： w(500)=5.1178%

在问题（二）中，m=0,从而rm*⊿t=0,其中v=4.871e12,r=3.6635132e10,w(t)=5%*wo 代入表示式（1）中，得：t=398.3(天)。

附录：

用MATLAB解问题（一）过程如下：

r=4.45365e10;v=5.941e12;wo=10;m=5;t=500; w=m+(wo-m).*exp(-r.*t) 结果： w =5.1178 问题（二）

r=3.6635132e10;v=4.871e12;m=0;w=0.05wo; w=m-(wo-m)*exp(-r*t)

第 32 页 共 46 页

a=log(w);y=a.*(-v)/r 结果：

y = 398.3120

七、插值与拟合

（一）实验目的

1.掌握插值与拟合的内容。

2.学会利用软件求解插值和拟合问题。

（二）实验课时

7.1在化工生产中常常需要知道丙烷在各种温度T和压力P下的导热系数K．下面是实验得到的一组数据：

T(℃) 68 68 87 87 106 106 140 140 p(103kN/m2) 9.7981 13.324 9.0078 13.355 9.7918 14.277 9.6563 12.463 K 0.0848 0.0897 0.0762 0.0807 0.0696 0.0753 0.0611 0.0651 试求T=99℃和P=10.3×103kN/m2下的K． p1=[9.7891,13.324]; k1=[0.0848,0.0897]; p2=[9.0078,13.355]; k2=[0.0762,0.0807]; p3=[9.7918,14.277]; k3=[0.0696,0.0753]; p4=[9.6563,12.463]; k4=[0.0611,0.0651]; a2=polyfit(p2,k2,1); a3=polyfit(p3,k3,1); x1=polyval(a2,10.3); x2=polyval(a3,10.3);

plot(10.3,x1,'k+',10.3,x2,'k+',p1,k1,p2,k2,p3,k3,p4,k4) xlabel('丙烷压力P(*10^3kPa)') ylabel('丙烷导热系数K')

title('在不同温度下丙烷导热系数与压力的关系图')

第 33 页 共 46 页

在T=87℃和T=106℃之间，仍采用线性近似来求T=99℃时的导热系数K。 x=[87 106]; y=[x1,x2];

z=polyval(a,99) a=polyfit(x,y,1); plot(99,z,'k+',x,y) grid

xlabel('丙烷温度T(℃)') ylabel('丙烷导热系数K')

title('压力10.3*10^3KPa时丙烷导热系数与温度的关系图')

所以，T=99℃和P=10.3*10^3kPa下的K=0.0729

7.2弹簧在力F的作用下伸长x，一定范围内服从胡克定律：F与x成正比，即F=kx．现在

第 34 页 共 46 页

得到下面一组F、x数据，并在(x，F)坐标下作图，可以看到当F大到一定数据值后，就不服从这个定律了．试由数据确定k，并给出不服从胡克定律时的近似公式． x 1 2 4 7 9 12 13 15 17 F 1.5 3.9 6.6 11.7 15.6 18.8 19.6 20.6 21.1

x=[0 1 2 4 7 9]；

f=[0 1.5 3.9 6.6 11.7 15.6]； A=polyfit(x，f，1) z=polyval(A，x)；

plot(x，f，'k+'，x，z，'r') 运行结果： A =

1.7085 0.0008

即F=1.7085x，说明当x≤9时，大致服从胡克定律．当x≥9后可以用如下二次函数来表示： 运行代码：

x=[9 12 13 15 17];

f=[15.6 18.8 19.6 20.6 21.1]; A=polyfit(x,f,2); z=polyval(A,x); plot(x,f,'k+',x,z,'b') 运行结果： A =

-0.0764 2.6728 -2.2613

第 35 页 共 46 页

八、概率统计、回归分析模型

（一）实验目的

1. 复习数据的录入、保存和调用 2. 直观了解统计基本内容.

3. 掌握用数学软件包Matlab求解概率、统计问题.

4. 综合提高用Matlab软件解决实际问题的能力. 学会分析问题、建立问题的数学表达式并加以求解和推广,体会数学建模的整个过程.

（二）实验课时

8.1 某校60名学生的一次考试成绩如下:

93 75 83 93 91 85 84 82 77 76 77 95 94 89 91 88 86 83 96 81 79 97 78 75 67 69 68 84 83 81 75 66 85 70 94 84 83 82 80 78 74 73 76 70 86 76 90 89 71 66 86 73 80 94 79 78 77 63 53 55 1)计算均值、标准差、方差、偏度、峰度,画出直方图; 2)检验分布的正态性;

3)若检验符合正态分布,估计正态分布的参数并检验参数计算均值、标准差、方差、偏度、峰度,画出直方图;

解： 1)

A=[93 75 83 93 91 85 84 82 77 76 77 95 94 89 91 88 86 83 96 81 79 97 78 75 67 69 68

84 83 81 75 66 85 70 94 84 83 82 80 78 74 73 76 70 86 76 90 89 71 66 86 73 80 94 79 78 77 63 53 55];

mean=mean(A); std=std(A); var=var(A);

skewness=skewness(A); kurtosis=kurtosis(A); hist(A,60)

disp([' mean = ',num2str(mean)]) disp([' std = ',num2str(std)]) disp([' var = ',num2str(var)])

disp([' skewness = ',num2str(skewness)]) disp([' kurtosis = ',num2str(kurtosis)])

mean = 80.1 std = 9.7106 var = 94.2949

skewness = -0.46818 kurtosis = 3.1529

第 36 页 共 46 页

2)检验分布的正态性; 解：

A=[93 75 83 93 91 85 84 82 77 76 77 95 94 89 91 88 86 83 96 81 79 97 78 75 67 69 68

84 83 81 75 66 85 70 94 84 83 82 80 78 74 73 76 70 86 76 90 89 71 66 86 73 80 94 79 78 77 63 53 55];

normplot(A)

学生的一次考试成绩近似服从正态分布.

3)若检验符合正态分布,估计正态分布的参数并检验参数. 解：

A=[93 75 83 93 91 85 84 82 77 76 77 95 94 89 91 88 86 83 96 81 79 97 78 75 67 69 68

84 83 81 75 66 85 70 94 84 83 82 80 78 74 73 76 70 86 76 90 89 71 66 86 73 80 94 79 78 77 63 53 55];

[muhat,sigmahat,muci,sigmaci]=normfit(A) [h,sig,ci]=ttest(A,80.1)

muhat =

80.1000

第 37 页 共 46 页

sigmahat =

9.7106

muci =

77.5915 82.6085

sigmaci =

8.2310 11.8436

估计出学生的一次考试成绩的均值为80.1,方差9.7106,均值的0.95置信区间为

[77.5915,82.6085],方差的0.95置信区间为[8.2310,11.8436].

h =

0 sig =

1 ci =

77.5915 82.6085

检验结果: 1. 布尔变量h=0, 表示不拒绝零假设,说明提出的假设学生的一次考试成

绩均值80.1是合理的.

2. 95%的置信区间为 [77.5915,82.6085], 它完全包括80.1, 且精度很高. 3. sig值为1, 远超过0.5, 不能拒绝原假设. 8.2 汽油价格

据说某地汽油的价格是每加仑115美分,为了验证这种说法,一位学者开车随机选择了一些加油站,得到某年一月和二月的数据如下：

一月：119 117 115 116 112 121 115 122 116 118 109 112 119 112 117 113 114 109 109 118 二月：118 119 115 122 118 121 120 122 128 116 120 123 121 119 117 119 128 126 118 125 1）分别用两个月的数据验证这种说法的可靠性; 2）分别给出1月和2月汽油价格的置信区间; 3）给出1月和2月汽油价格差的置信区间;

>> a=[119 117 115 116 112 121 115 122 116 118 109 112 119 112 117 113 114 109 109

118]

a = Columns 1 through 14

119 117 115 116 112 121 115 122 116 118 109 112 119 112 Columns 15 through 20

117 113 114 109 109 118

>> b=[118 119 115 122 118 121 120 122 128 116 120 123 121 119 117 119 128 126 118

125]

b =

Columns 1 through 14

118 119 115 122 118 121 120 122 128 116 120 123 121 119 Columns 15 through 20

117 119 128 126 118 125 >> [n,x]=HIST(a)

第 38 页 共 46 页

Warning: Function call HIST invokes inexact match

C:\\MATLAB7\\toolbox\\matlab\\datafun\\hist.m.

n =

3 0 3 2 2 2 4 2 0 2 x =

Columns 1 through 8

109.6500 110.9500 112.2500 113.5500 114.8500 116.1500 117.4500

118.7500

Columns 9 through 10 120.0500 121.3500 >> [n,x]=HIST(b) n =

2 1 3 5 2 2 1 1 1 2 x =

Columns 1 through 8

115.6500 116.9500 118.2500 119.5500 120.8500 122.1500 123.4500

124.7500

Columns 9 through 10 126.0500 127.3500 >> [nuhat,muci]=expfit(a) nuhat =

115.1500 muci =

77.6183 188.5152 >> [nuhat, muci]=expfit(b) nuhat =

120.7500 muci =

81.3930 197.6831

>> [muhat,sigmahat,muci,sigmaci]=normfit(a) muhat = 115.1500 sigmahat = 3.8699 muci = 113.3388 116.9612 sigmaci = 2.9430 5.6523

>> %muhat u的估值muci为参数u的置信区间估计的，sigmahat 为6的估计值，

sigmaci 6的置信区间

>> [muhat,sigmahat,muci,sigmaci]=normfit(b) muhat = 120.7500 sigmahat = 3.7116 muci = 119.0129 122.4871 sigmaci = 2.8227 5.4211

8.3根据牙膏销售量与销售差价、广告费用等数据（见下表）,建立数学模型,分析牙膏销售与其它因素的关系,为制订价格策略和广告投入策略提供数量依据.

第 39 页 共 46 页

销售周期 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 广告费用（百万元） 5.5 6.75 7.25 5.5 7 6.5 6.75 5.25 5.25 6 6.5 6.25 7 6.9 6.80 6.8 7.1 7 6.8 6.5 6.25 6 6.5 7 6.8 6.8 6.5 5.75 5.8 价格差（元） -0.05 0.25 0.60 0 0.25 0.2 0.15 0.05 -0.15 0.15 0.2 0.1 0.4 0.45 0.35 0.3 0.5 0.5 0.4 -0.05 -0.05 -0.1 0.2 0.1 0.5 0.6 -0.05 0 0.05 销量量（百万支） 7.38 8.51 9.52 7.5 9.33 8.28 8.75 7.87 7.1 8 7.89 8.15 9.1 8.86 8.9 8.87 9.26 9 8.75 7.95 7.65 7.27 8 8.5 8.75 9.21 8.27 7.67 7.93 30 6.8 0.55 9.26 解：

模型假设

假设该公司的价格差与广告费用是相互独立的，互不影响。

设：销售量为y，价格差x1，广告费x2

y=[7.38 8.51 9.52 7.5 9.33 8.28 8.75 7.87 7.1 8 7.89 8.15 9.1 8.86 8.9 8.87 9.26 9 8.75

7.95 7.65 7.27 8 8.5 8.75 9.21 8.27 7.67 7.93 9.26];

y=y.';

x_1=[-0.05 0.25 0.60 0 0.25 0.2 0.15 0.05 -0.15 0.15 0.2 0.1 0.4 0.45 0.35 0.3 0.5 0.5 0.4

-0.05 -0.05 -0.1 0.2 0.1 0.5 0.6 -0.05 0 0.05 0.55];

x_1=x_1.';

x_2=[5.5 6.75 7.25 5.5 7 6.5 6.75 5.25 5.25 6 6.5 6.25 7 6.9 6.80 6.8 7.1 7 6.8 6.5 6.25 6

6.5 7 6.8 6.8 6.5 5.75 5.8 6.8];

x_2=x_2.'; x_4=(x_2).^2;

第 40 页 共 46 页

eye=ones(1,30); eye=eye.';

x=[eye x_1 x_2 x_4];

[b,bint,r,rint,stats]=regress(y,x,0.05) b =

17.3244 1.3070 -3.6956 0.3486

bint =

5.7282 0.6829 -7.4989 0.0379 r =

-0.0988 -0.0795 -0.1195 -0.0441 0.4660 -0.0133 0.2912 0.2735 -0.2351 0.1031 -0.4033 0.1747 0.0400 -0.1504 0.1284 0.1637 -0.0527 -0.1907 -0.0870 -0.0165 -0.1292 -0.3002 -0.2933 -0.1679 -0.2177 0.1116 0.3035 0.0693 0.2474 0.2270 rint =

-0.5270 -0.5309 -0.5106 -0.4731 0.0813 28.9206 1.9311 0.1077 0.6594 0.3294 0.3718 0.2716 0.3848 0.8507

第 41 页 共 46 页

-0.4609 0.4343 -0.1374 0.7197 -0.0870 0.6340 -0.5960 0.1258 -0.3280 0.5341 -0.8190 0.0125 -0.2618 0.6112 -0.4032 0.4832 -0.5933 0.2925 -0.3207 0.5775 -0.2841 0.6116 -0.4830 0.3776 -0.6248 0.2434 -0.5348 0.3609 -0.4423 0.4092 -0.5609 0.3024 -0.7181 0.1177 -0.7243 0.1377 -0.5548 0.2190 -0.6449 0.2095 -0.2994 0.5226 -0.1037 0.7106 -0.3714 0.5099 -0.1807 0.6755 -0.1890 0.6430

stats =

0.9054 82.9409 0.0000 0.0490 参数 参数估计值 // 17.3244 // 1.3070 // -3.6956 // 0.3486 R//=0.9054 F=82.9409 p<0.0001 s//=0.0490 y=17.3244+1.3070x1-3.6956x2+0.3486

参数置信区间 [5.7282，28.9206] [0.6829，1.9311] [-7.4989，0.1077 ] [ 0.0379，0.6594] 九、计算机模拟

（一）实验目的

1.掌握计算机模拟的内容。

2.学会利用软件求解计算机模拟问题。

（二）实验课时

9.1.当天生产的产品必须售出，否则就会变质。该产品单位成本为 2.5 元，单位产品售价为5元。企业为避免存货过多而造成损失，拟从以下两种库存方案中选出一个较优的方案：

第 42 页 共 46 页

方案甲：按前1天的销售量作为当天的库存量；

方案乙：按前2天的平均销售量作为当天的库存量。 假定市场对该产品的每天需求量是一个随机变量，但从以往的统计分析得知它服从正态分布N(135 , 22.42)。

function [TL1,TL2]=kucum(T,S1,S21,S22)TL1=0;TL2=0;N=1;while N<TQ1=S1;Q2=(S21+S22)/2;D=normrnd(135,22.4);if D<Q1; S3=Q1;else S3=D;endif D<Q2 S4=Q2;elseS4=D;endL1=5*S3-2.5*Q1;L2=5*S4-2.5*Q2;TL1=TL1+L1;TL2=TL2+L2;N=N+1;S1

=S3;S22=S21;S21=S4;end TL1<

TL2，所以方案乙较优

9.2.某报童以每份0.03元的价格买进报纸，以0.05元的价格出售．根据长期统计，报纸每天的销售量及百分率为

200 210 220 230 240 250 销售量

0.10 0.20 0.40 0.15 0.10 0.05 百分率

第 43 页 共 46 页

已知当天销售不出去的报纸，将以每份0.02元的价格退还报社．试用模拟方法确定报童每天买进报纸数量，使报童的平均总收入为最大？

主文件main.m：

BUYMIN=200; % 每天的最小购买量 BUYMAX=250; % 每天的最大购买量

SIMUDAY=1.0e+5; % 模拟时间 sell_amount=200:10:250; % 销售量 percentage=[0.1 0.3 0.7 0.85 0.95 1]; % 百分率 buy_amount=0; ave_profit=0;

for loop_buy=BUYMIN:BUYMAX sum_profit=0;

for loop_day=1:SIMUDAY

index=find(percentage>=rand); % 产生随机数，用于决定当天的销售量 sum_profit=sum_profit+GetProfit(loop_buy,sell_amount(index(1))); end

buy_amount=[buy_amount,loop_buy]; % 循环嵌套 ave_profit=[ave_profit,sum_profit/SIMUDAY]; % 循环嵌套 end

buy_amount(1)=[]; % 第一个元素置空 ave_profit(1)=[];

[val,id]=max(ave_profit) % 显示最大平均收入val

buy=buy_amount(id) % 显示在平均收入最大情况下的每天的购买量buy xlabel='每天的购买量'; ylabel='平均利润';

plot(buy_amount,ave_profit,'*:'); % 函数GetProfit.m代码： function re=GetProfit(a,b)

if a

else % 供过于求：报童购买量大于销售量 re=b*(0.05-0.03)+(a-b)*(0.02-0.03); end 运行结果：

第 44 页 共 46 页

val =4.2801 id =21 buy = 220

该结果说明当报童每天买进报纸数量为220，报童的平均总收入为最大，且最大为

4.2801。

第 45 页 共 46 页

第 46 页 共 46 页