数学模型实验

数学模型实验
实验一
一：商人过河
1、问题提出：三名商人各带一个随从乘船渡河，一只小船只能容纳二人，由他们自己划行。

随从们密约，在河的任一岸，一旦随从人数比商人多，就商人越货。

但是如何乘船渡河的大权掌握在商人们中。

商人们怎样安全渡河呢？
2、模型构成 记第k 次渡河前此案的商人数为k x ，随从数
.3,2,1,0,,,2,1,==k k k y x k y 将二维向量),(k k k y x s =定义为状态。

安全渡河条件下的状态集合称为允许状态集合，记作S .
}2,1;3,2,1,0,3;,2,1,0,0|),{(=======y x y x y x y x S (1)不难验证，S 对此岸和彼岸都是安全的.
记第k 次渡船上的商人数为k u ，随从数为k v .将二维向量),(k k k v u d =定义为决策，允许决策集合记作 ，由小船的容量可知
}2,1,0,.21|),{(=≤+≤=v u v u v u D （2） 因为k 为奇数时船从此岸驶向彼岸，k 为偶数时船由彼岸驶回此岸，所以状态随决策k d 变化的规律是
k k k k d s s )1(1-+=+ （3）
（3）式称为状态转移律.制定安全渡河方案归结为 多步决策模型：
求决策),,2,1(n k D d k =∈，使状态S s k ∈按照转移规律（3），由初始状态)3,3(1=s 经有限步n 到达状态)0,0(1=+n s
3、程序设计：
function jueche=guohe
n=input('输入商人数目: ');
nn=input('输入仆人数目: ');
nnn=input('输入船的最大容量: ');
if nn>n
n=input('输入商人数目:');
nn=input('输入仆人数目:');
nnn=input('输入船的最大容量:');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 决策生成
jc=1; %决策向量放在矩阵d 中，jc 为插入新元素的行标初始为1；
for i=0:nnn
for j=0:nnn
if (i+j<=nnn)&(i+j>0) % 满足条D={(u,v)|1<=u+v<=nnn,u,v=0,1,2} d(jc,1:3)=[i,j,1];%生成一个决策向量立刻扩充为三维；
d(jc+1,1:3)=[-i,-j,-1]; % 同时生成他的负向量；
jc=jc+2; % 由于生成两个决策向量，则jc要向下移动两个；
end
end
j=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 状态数组生成
kx=1; % 状态向量放在A矩阵中，生成方法同矩阵生成；
for i=n:-1:0
for j=nn:-1:0
if ((i>=j)&((n-i)>=(nn-j)))|((i==0)|(i==n))
% (i>=j)&((n-i)>=(nn-j)))|((i==0)|(i==n))为可以存在的状态的约束条件
A(kx,1:3)=[i,j,1]; %生成状态数组集合D `
A(kx+1,1:3)=[i,j,0];
kx=kx+2;
end
end
j=nn;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 将状态向量生
成抽象矩阵
k=(1/2)*size(A,1);
CX=zeros(2*k,2*k);
a=size(d,1);
for i=1:2*k
for j=1:a
c=A(i,:)+d(j,:) ;
x=find((A(:,1)==c(1))&(A(:,2)==c(2))&(A(:,3)==c(3))) ;
v(i,x)=1; %x为空不会改变v值
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dijstra算法x=1; y=size(A,1);
m=size(v,1);
T=zeros(m,1);
T=T.^-1;
lmd=T;
P=T;
S=zeros(m,1);
S(x)=1;
P(x)=0; lmd(x)=0;
k=x;
while(1)
a=find(S==0);
aa=find(S==1);
if size(aa,1)==m
break;
end
for j=1:size(a,1)
pp=a(j,1);
if v(k,pp)~=0
if T(pp)>(P(k)+v(k,pp))
T(pp)=(P(k)+v(k,pp));
lmd(pp)=k;
end
end
end
mi=min(T(a));
if mi==inf
break;
else
d=find(T==mi);
d=d(1);
P(d)=mi;
T(d)=inf;
k=d;
S(d)=1;
end
end
if lmd(y)==inf
jueche='can not reach';
return;
end
jueche(1)=y;
g=2; h=y;
while(1)
if h==x
break;
end
jueche(g)=lmd(h);
g=g+1;
h=lmd(h);
end
jueche=A(jueche,:);
jueche(:,3)=[];
在matlab窗口输入
商人数目：3
随从数目：3
船的最大容量：2
二、施救药物中毒
1、问题提出：一个孩子两个 小时前一次性吞了11片治疗哮喘病的、剂量为没片100mg 的氨茶碱片，已经出现呕吐、头晕等不良症状.氨茶碱的每次用量成人是100-200mg/kg ，儿童是3-5mg/kg ，如果用量过高，会使血液浓度过高。

当血药浓度达到100ug/ml 时，会出现严重中毒，达到200则可致命.
2、问题调查与分析：人体服用一定量的药物后，血药浓度与人体的血液总量有关。

血液总量约为体重的7%-8%，即体重50-60kg 的成年人有4000ml 左右的血液。

孩子的体重约为成年人的一半，所以血液总量为2000ml 。

血液系统对药物吸收和排除率的半衰期分别为5h 、6h 。

3、模型的假设和建立：为了判断孩子的血液浓度会不会达到危险的水平，需要求得胃肠道和血液系统中的药量随时间变化的规律。

记肠胃道的药量为x(t),血液系统的药量为y(t)，时间t 以孩子误服药的时刻为起点（t=0）。

根据前面的调查与分析可作以下假设：
（1）胃肠道中药物向血液系统的转移率与药量X （t ）成正比，比例系数记作 （>0），总剂量x mg 的药物在t=0瞬间进入胃肠道。

（2）血液系统中药物的排除率与药量y （t ）成正比，比例系数记作 （u>0），t=0是血液中无药物。

（3）氨茶碱被吸收的半衰期为5小时，排除的半衰期6小时。

孩子的血液总量为2000ml ，成人的血液总量为4000ml 。

根据假设对胃肠道中药量X （t ）和血液系统中药量y （t ）建立服下模型。

有假设1，X （0）=x ，随着药物从胃肠道向血液系统的转移，X （t ）下降的速度与X （t ）本身成正比（比例系数 0>λ）,所以X （t ）满足微分方程
1100)0(,=-=x x dt
dx λ （1） 由假设2，y （0）=0，药物从胃肠道向血液系统的转移相当于血液系统对药物的吸收，y （t ）由于吸收作用而增长的速度是x ，由于排除而减少的速度与y
（t ）本身成正比（比例系数 >0），所以y （t ）满足微分方程
0)0(,=-=y y x dt
dy μλ （2） 方程（1），（2）中的参数 和可由假设3的半衰期确定。

模型求解 微分方程（1）是可分离变量方程，容易得到
t e t x λ-=1100)( （3）
表明胃肠道中的药量X （t ）随时间单调减少并趋于零。

为了确定，利用药
物吸收的半衰期为5，即2/11002/)0(1100)(5=-==-x e t x t ，得)/1(1386.05/)2(l n h ==λ。

因此可求得)(1100)(t t e e t y λμμ
λλ----= （4）将1386.0=λ和1155.0=μ带入（3）、（4）中可得到想
t e t x 1386.01100)(-= （5）
)(1100)(1386.01155.0t t e e t y ----=μ
λλ （6） 根据假设4，孩子的血液总量为2000ml ，出现严重中毒的血药浓度100mg g /μ和致命的血药浓度200mg g /μ分别相当与血液中药量y 达到200mg 和400mg 。

由（6）容易精确的算出y （2）=236.5mg ，计算药量达到400mg 的时间为t1，解非线性方程
400)(6600116.138.01155.0=---t t e e ，用matlab 软件计算可以得到h t 87.41=.
施救方案：若用口服活性炭来吸附药物的办法施救，药物的排除率可增加到μ的2倍，即0.2310.设孩子达到医院时刻（t=2）就开始施救，由（2）、（3），建立新的模型为
5.236)2(,1100,2,==≥-=-z x t z x dt
dy t λμλ （7） 解得2310.0-=μ时，（7）的解为
2,5.16091650)(2310.01386.0≥-=--t e e t z t t （8）
用matlab 对x(t),y(t ），z （t ）作图得：
施救后血液系统中药量)(t z 。

)(),(t y t x 。

由图中可以看出血液中药量z （t ）达到最大值的时间约在t=5h ，即到达医院施救后3h ，其精确值可由（8）算出，记t3，t3=5.26h ，且z(t3)=318.4mg ，远低于y(t)的最大值和致命水平。

4、程序如下：
t=0:25;
x=1100*exp(-0.1386*t);
plot(t,x);
hold on
y=6600*(exp(-0.1155*t)-exp(-0.1386*t));
plot(t,y)
hold on
t=2:25
z=1650*exp(-0.1386*t)-1609.5*exp(-0.2310*t);
plot(t,z)
title('胃肠道中药量x （t ）和血液系统中药量y （t ）')
xlabel('t/h')
ylabel('x,y,z/mg')
text(5,400,'y(t)')
text(2.5,800,'x(t)')
text(5,300,'z(t)')
grid on
三、饮料厂的生产与检修
（一）饮料厂的检修与生产计划
1、问题提出：如果工厂必须在未来四周的某一周中安排一次设备检修，检修将占用当周的15千箱的生产能力，但会使检修以后每周的生产能力提高5千箱，则检修应该安排在哪一周？
2、模型假设：设饮料厂在第一周开始没有库存；从费用最小考虑，自然的假设第4周末也不能有库存；周末有库存时需支出一周的存储费；每周末的库存量就是下周初的库存量。

3、模型构造与求解：
列出线性规划模型与约束条件：
)(2.05.54.51.50.5min 3214321y y y x x x x z ++++++= (1)
S.t 1511=-y x (2-1)
25211=-+y y x (2-2) 35323=-+y y x (2-3) 2534=+y x (2-4) 20,45,40,304321≤≤≤≤x x x x (3) 0,,,,,,3214321≥y y y x x x x (4)
将以上模型输入LINGO 求解，可以得到最优解)5,15,0,20,25,40,15(),,,,,,(3214321=y y y x x x x ,
实际上，由于检修后生产能力发生变化，把检修安排在第一周和第二周不一定是最佳选择。

引入0-1变量4321,,,ωωωω，用1=t ω表示检修安排在t 周)4,3,2,1(=t 。

所以上面的约束条件需修改为
301511≤+ωx （1）
401515122≤-+ωωx （2）
4555151233≤--+ωωωx （3）
205551532144≤---+ωωωωx （4）
同时，模型中还需要增加关于4321,,,ωωωω的约束
14321=+++ωωωω （5）
}1,0{,,,4321∈ωωωω （6）
利用ＬＩＮＧＯ编辑程序及结果如下：
线性规划：
程序（１）
结果：
程序（2）
即当将检修安排在第一周，每周的生产量分别为15千箱、45千箱、15千箱、25千箱，最小总费用从528千元下降为527千元。

（二）饮料厂的生产批量问题
1、问题提出：某饮料厂使用同一条生产线轮流生产多种饮料满足市场需求。

如果某周开工生产其中一种饮料，就要清洗设备和更换部分部件，于是需支出生产准备费8千元。

假设其未来四周需求量、生产能力、生产成本与例1给出的完全相同。

问应如何安排这种饮料的生产计划，在按时满足市场需求的条件下，使生产该饮料的总费用最小？
2、问题分析假设：这一问题与例1中最主要的差别是：除了考虑随产品数量的变化的费用外，还要考虑与产品数量无关的费用，即生产准备费。

3、设需要考虑时间跨度T 个小时段，在t 时段的市场需求dt ，生产能力为Mt （t=1,2,...T ）。

如果第t 时段开工生产，则需付出生产准备费为0≥t s ，单件生产成本为0≥t c 。

在t 时段末，如果有产品库存，单件产品1个时段的存储费为0≥t h 。

假设在t 时段，产品的生产量为)0(≥t x ,t 时段末产品的库存为)0(≥t y ，合理的假设00==τy y 。

产量、需求与库存之间的平衡关系同前。

为了表述在t 时段是否生产这种饮料，从而确定是否要支付生产准备费，引进0-1变量t ω，
1=t ω表示生产，0=t ω表示不生产。

所以模型如下：
∑=++=T
t t t t t t t y h x c s z 1)(min ω （7）
..t s ,1t t t t d y x y =-+- T t ,,2,1 = （8）
T t x x t t t ,,2,1,0,0,0,1{ ==>=ω （9）
T t M x t t ,,2,1, =≤ （10） 00==τy y ， （11） T t y x t t ,,2,1,0, =≥ （12）
将（9），（10）合并表示为t t t M x ω≤，或用T t M x t t t ,,2,1,0 =≤-ω代替。

用LINGO 写程序求解如下：
程序（３）
求得生产计划为：4周的生产量分别为15、40、45、0千箱，最小总费用为554千元。

与例1比，这里只在前三周生产，虽然存储费增加了，但节省了生产准备费。

四、钢管和易拉罐下料
（一）钢管下料
1、问题提出：某钢管零售商从钢管厂进货。

按照顾客的要求切割后出售，从钢管厂进货时得到的原料都是19m 。

（1）现有一个客户需要50根4m 、20根6m 和15根8m 的钢管。

应如何下料最节省？ （2）零售商如果采用不同切割模式太多，将会导致生产过程的复杂化，从而增加生产和管理成本，所以该零售商规定采用的不同模式不能超过3种。

此外，该客户除了需要（1）中的三种钢管，还需要10根5m 的钢管。

应如何下料最节省？
2、问题分析与假设：决策变量 用x 表示按照第i 种模式)7,2,1( =i 切割的原料钢管的根数，显然它们应该是非负数。

3、模型建立： 决策目标与约束条件：
765432113333min x x x x x x x z ++++++= (1)
76543212min x x x x x x x z ++++++= (2)
5023454321≥++++x x x x x (3) 20326542≥+++x x x x (4) 152753≥++x x x (5)
对于问题二，与一类似，一个合理的切割模式的余料不应该大于或等于客户需要的钢管的最小尺寸，切割计划只适合合理的切割模式，而由于本题中的参数都是整数，所以合理的切割模式余量不能大于3.所以可以用i x 表示第i 中模式（i=1,2,3）切割的原料钢管的根数，显然它们应当是非负整数，设所使用的第i 中切割模式下没根原料钢管生产4m,5m,6m 和8m 的钢管数量分别为i i i i r r r r 4321,,,（非负整数）。

决策目标和约束条件：
321min x x x z ++= （6） 50313212111≥++x r x r x r （7） 10323222121≥++x r x r x r （8） 20333232131≥++x r x r x r （9） 153********≥++x r x r x r （10）
每一种切割模式必须可行、合理，所以每根原料钢管的成品量不能超过19m ，也不能少于16m （余量不能大于3m ），于是
1986541641311211≤+++≤r r r r （11）
1986541642322212≤+++≤r r r r （12） 1986541643332313≤+++≤r r r r （13）
由于三种切割模式的排列顺序是无关紧要的，所以不妨增加一些约束条件。

321x x x ≥≥ （14）
原料钢管的总根数不可能少于26119158206105504=+⎥⎦
⎤
⎢⎣⎡
⨯+⨯+⨯+⨯根。

其次，考虑一种非常特殊的生产计划：第一种切割模式下只生产4m 钢管，一根原料钢管切割成4根4m 钢管，为满足50根4m 钢管的需求，需要13根原料钢管;第二种切割模式下只生产5m 、6m 钢管，一根原料钢管切割成1根5m 钢管和2根6m 钢管，为满足10根5m 和20根6m 钢
管需求，需要10根原料钢管；第三中切割模式下只生产8m 钢管，一根原料钢管切割成2根8m 钢管，为满足15根8m 钢管的需求，需要8根原料钢管。

于是满足要求的这种生产计划共需3181013=++根原料钢管，这就得到了最优解的一个上界。

所以可增加以下约束条件：
3126321≤++≤x x x （15）
4、程序设计：
程序（1）：
结果：
解得最优解为：15,1252==x x ，其余全为0,。

即按照模式2切割12根，模式5切割15根，共27根，总余料量为27根。

显然在总余料量最小的目标下，最优解将是使用余料尽可能小的切割模式，这将会导致切割原料钢管的总根数较多。

程序（2）：
结果：
即按照模式1、2、3分别切割10,10,8根原料钢管，使用原料总数为28根。

第一种切割模式下一根原料钢管切割成3根4m 钢管和1根6m 钢管；第二中切割模式下一根原料钢管切割成2根4m 、1根5m 钢管和1根6m 钢管；第三种切割下一根原料钢管切割成2根8m 钢管。

但这个模型的解并不唯一。

（二）易拉罐下料
1、问题提出：某公司采用一套冲压设备生产一种罐装饮料的易拉罐，易拉罐为圆柱形，包括罐身、上盖和下底，罐身高10cm ，上盖和下底的直径均为5cm 。

该公司使用两种不同规格的镀锡板原料：规格1的镀锡板为正方形，边长24cm ；规格2的镀锡板为长方形，长、宽分别为32cm 和28cm 。

由于生产设备和生产工艺的限制，对于规格2的镀锡板原料。

只能按照模式4进行冲压。

使用模式1、模式
2、模式
3、模式4 进行每次冲压所需要的时间分别为1.5s 、2s 、1s 、3s 。

该工厂每周工作40h ，每周可供使用的规格1、规格2的镀锡板原料分别为5万张和2万张。

目前每只易拉罐的利润为0.10元，原料余料损失为0.001元/cm(如果周末有罐身、上盖或下底不能配套组装成易拉罐出售，也看做是原料余料损失)。

2、问题分析；已知上盖和下底的直径d=5cm ，可得面积为2
26.19/cm d d s ≈=π，周长为
cm d L 7.15≈=π；已知罐身高为h=10cm 。

可得其面积为21.157cm hL S ≈=。

于是模式1
下的余料损失为226.2221024cm S s ≈--。

3、模型建立：
用i x 表示按照第i 种模式的冲压次数（i=1,2,3,4），1y 表示一周生产的易拉罐个数。

为计算不能配套组装的罐身和底、盖造成的原料损失，用2y 表示不配套的罐身个数，3y 表示不配套的底、盖个数。

虽然实际上i x 和321,y y y 应该是整数。

但是由于生产量相当大，可以把它们看成是实数，从而用线性规划模型处理。

决策目标和约束条件：
假设每周生产的易拉罐能够全部售出，公司每周的销售量利润是0.1y1.原料余料损失包括两部分：4种冲压模式下的余料损失和不配套的罐身和底、盖造成的原料损失。

按照前面的计算及前面的结果，总损失为
().6.191.1575.1698.2613.1836.222001
.0324321y y x x x x +++++ ()32432116.191.1575.1698.2613.1836.222001
.01.0max y y x x x x y z +++++-= （16）
144000325.14321≤+++x x x x （17）
每周可供使用的规格1、规格2的镀锡板原料分别为50000张和20000张，即
50000321≤++x x x （18）
200004≤x （19）
一周生产的罐身个数为42142x x x ++，一周生产的底、盖个数为4321516410x x x x +++，因为应尽可能将它们配套组装成易拉罐销售。

所以1y 满足
}2/)516410(,42min{4321421x x x x x x x y +++++= （20）
这时不配套的罐身个数2y 和不配套的底、盖个数2y 应为
1421242y x x x y -++= （21） 1432132516410y x x x x y -+++= （22）
421142x x x y ++≤ （23） 2/)516410(43211x x x x y +++≤ （24）
由于（17）——（19）中右端的数值过大，模型中数据之间的数量级不匹配，此时ＬＩＮＧＯ在计算时容易产生较大的误差。

因此把它改为：
4.1442
5.14321≤+++x x x x （25） 5321≤++x x x （26）
24≤x （27）
4、程序设计：
实验二
一、线性回归
1、问题提出：在某产品表明腐蚀刻线，下表是试验活得的腐蚀时间（x ）与腐蚀深度（y ）间的一组数据。

试研究两变量（x ，y ）之间的关系。

X 5 5 10 20 30 40 50 60 65 90 100 y
5
6
8
13
15
17
19
25
25
29
35
其中：x------腐蚀时间(秒)； ------腐蚀深度（y ）(m μ)。

要求： 1）画出散点图，并观察y 与x 的关系；
2）求y 关于x 的线性回归方程： y a
bx =+ ，求出a 与b 的值； 3）对模型和回归系数进行检验；
4）预测x=120时的y 的置信水平为0.95的预测区间。

5）编程实现上述求解过程。

程序解释：
[b, bint,r,rint,stats]=regress(Y,X,alpha) 求回归系数的点估计和区间估计、并检验回归模型
b ：回归系数的估计值
bint ：表示回归系数的区间估计. r ：表示残差
rint ：表示置信区间
stats ：表示用于检验回归模型的统计量,有三个数值：相关系数r2、F 值、与F 对应的概率p
2、程序设计：
>> x1=[0 5 10 20 30 40 50 60 65 90 100]'; >> x=[ones(11,1),x1];
>> y=[5 6 8 13 15 17 19 25 25 29 35]'; >> [b,bint,r,rint,stats]=regress(y,x,0.05) b =
5.6273 0.2874
bint =
4.0401 7.2146 0.2578 0.3171 r =
-0.6273
-1.0646
-0.5018
1.6238
0.7493
-0.1251
-0.9996
2.1259
0.6887
-2.4974
0.6281
rint =
-3.4880 2.2333
-3.9122 1.7830
-3.5054 2.5018
-1.2128 4.4603
-2.3692 3.8678
-3.3232 3.0729
-4.0904 2.0912
-0.5307 4.7826
-2.3853 3.7627
-4.5073 -0.4875
-1.9681 3.2244
stats =
0.9816 479.8521 0.0000 1.9575 画图：>> rcoplot(r,rint)
结果分析：从残差图可以看出，除第10个数据外,其余数据的残差离零点均较近，且残差的置信区间均包含零点，第10个数据可视为异常点（剔除）。

所以剔除数据再计算：
clear all;clc
x1=[0 5 10 20 30 40 50 60 65 100]';
x=[ones(10,1), x1];
y=[5 6 8 13 15 17 19 25 25 35]';
[b,bint,r,rint,stats]=regress(y,x,0.05)
rcoplot(r,rint)
运行结果为：
b =
5.3232
0.3020
bint =
4.0806 6.5658
0.2763 0.3277
r =
-0.3232
-0.8333
-0.3434
1.6364
0.6162
-0.4040
-1.4242
1.5556
0.0455
-0.5253
rint =
-2.5344 1.8880
-3.0035 1.3368
-2.6625 1.9756
-0.3081 3.5809
-1.7762 3.0085
-2.8398 2.0317
-3.5243 0.6758
-0.4081 3.5192
-2.3014 2.3923
-2.2415 1.1910
stats =
0.9892 733.5467 0.0000 1.1080
从残差图可以看出，数据的残差离零点均较近，且残差的置信区间均包含零点，说明回归模型
y= 5.3232+0.3020x能较好的符合原始数据，
二、多元线性回归问题
1、问题提出：根据下述某猪场25头育肥猪4个胴体性状的数据资料，试进行瘦肉量y
对眼肌面积（x1）、腿肉量(x2)、腰肉量(x3)的多元线性回归分析。

序号瘦肉量
y(kg)
眼肌面积
x1(cm2)
腿肉量
x2(kg)
腰肉量
x3(kg)
序
号
瘦肉量
y(kg)
眼肌面积
x1(c m2)
腿肉量
x2(kg)
腰肉量
x3(kg)
1 15.0
2 23.7
3 5.49 1.21 1
4 15.94 23.52 5.18 1.98 2 12.62 22.34 4.32 1.3
5 15 14.33 21.8
6 4.86 1.59 3 14.86 28.84 5.04 1.92 16 15.11 28.95 5.18 1.3
7 4 13.9
8 27.67 4.72 1.4
9 17 13.81 24.53 4.88 1.39 5 15.91 20.83 5.35 1.56 18 15.58 27.65 5.02 1.66 6 12.47 22.27 4.27 1.50 19 15.85 27.29 5.55 1.70 7 15.80 27.57 5.25 1.85 20 15.28 29.07 5.26 1.82 8 14.32 28.01 4.62 1.51 21 16.40 32.47 5.18 1.75 9 13.76 24.79 4.42 1.46 22 15.02 29.65 5.08 1.70 10 15.18 28.96 5.30 1.66 23 15.73 22.11 4.90 1.81 11 14.20 25.77 4.87 1.64 24 14.75 22.43 4.65 1.82 12 17.07 23.17 5.80 1.90 25 14.35 20.04 5.08 1.53 13
15.40
28.57
5.22
1.66
要求：
1）画出散点图y 与1x ，y 与2x ，y 与3x 并观察y 与321,,x x x 的关系；
2）求y 关于321,,x x x 的线性回归方程： 0112233y a a x a x a x =+++-----（1），求出0123,,,a a a a 的值；
3）对上述回归模型和回归系数进行检验；
4）再分别求y 关于单个变量321,,x x x 的线性回归方程： 10111y a a x =+----（2）， 20222y a a x =+-----（3）, 30333
y a a x =+--- --（4）求出ij a 的值； 分别求y 关于两个变量x1，x2, x3的线性回归方程： 10111122y a a x a x =++----(2’)， 20211222y a a x a x =++---（3’）, 30311322
y a a x a x =++ --- --（4’）求出系数ij a 的值；并说明这六个回归方程对原来问题求解的优劣。

5）编程实现上述求解过程。

2、程序设计：
>> x1=[23.73 22.34 22.34 27.67 20.83 22.27 27.57 28.01 24.79 28.96 25.77 23.17 28.57 23.52 21.86 28.95 24.53 27.65 27.29 29.07 32.47 29.65 22.11 22.43 20.04];
>> x2=[5.49 4.32 5.04 4.72 5.35 4.27 5.25 4.62 4.42 5.30 4.87 5.80 5.22 5.18 4.86 5.18 4.88 5.02 5.55 5.26 5.18 5.08 4.90 4.65 5.08];
>> x3=[1.21 1.35 1.92 1.49 1.56 1.50 1.85 1.51 1.46 1.66 1.64 1.90 1.66 1.98 1.59 1.37 1.39 1.66 1.70 1.82 1.75 1.70 1.81 1.82 1.53];
>> y=[15.02 12.62 14.86 13.98 15.91 12.47 15.80 14.32 13.76 15.18 14.20 17.07 15.40 15.94 14.33 15.11 13.81 15.58 15.85 15.28 16.40 15.02 15.73 14.75 14.35]; >> subplot(2,2,1);plot(x1,y ,'b*');title('眼肌面积') >> subplot(2,2,2);plot(x2,y ,'r*');title('腿肉量')
>> subplot(2,2,3);plot(x3,y,'k*');title('腰肉量')
x1=[23.73 22.34 22.34 27.67 20.83 22.27 27.57 28.01 24.79 28.96 25.77 23.17 28.57 23.52 21.86 28.95 24.53 27.65 27.29 29.07 32.47 29.65 22.11 22.43 20.04]';
>> x=[ones(25,1),x1];
>> y=[15.02 12.62 14.86 13.98 15.91 12.47 15.80 14.32 13.76 15.18 14.20 17.07 15.40 15.94 14.33 15.11 13.81 15.58 15.85 15.28 16.40 15.02 15.73 14.75 14.35]';
>> [b,bint,r,rint,stats]=regress(y,x,0.05)
结果为：
b =
12.5432
0.0931
bint =
9.0744 16.0119
-0.0423 0.2284
r =
0.2680
-2.0026
0.2374
-1.1387
1.4280
-2.1461
0.6906
-0.8303
-1.0906
-0.0588
-0.7418
2.3702
0.1975
1.2076
-0.2479
-0.1278
-1.0164
0.4632
0.7667
0.0310
0.8345
-0.2830
1.1288
0.1190
-0.0585
rint =
-1.9234 2.4595
-3.9798 -0.0254
-1.9255 2.4004
-3.2653 0.9879
-0.5896 3.4456
-4.0916 -0.2005
-1.4751 2.8563
-2.9772 1.3166
-3.2429 1.0617
-2.2108 2.0933
-2.9240 1.4404
0.4515 4.2889
-1.9645 2.3595
-0.9180 3.3332
-2.3964 1.9006
-2.2796 2.0239
-3.1740 1.1411
-1.7128 2.6391
-1.3990 2.9324
-2.1176 2.1796
-1.1110 2.7801
-2.4076 1.8416
-0.9719 3.2296
-2.0482 2.2863
-2.1359 2.0189
stats =
0.0809 2.0244 0.1682 1.1342
画图：rcoplot(r,rint)
同理题一可知，剔除残差较大的数据。

>> x1=[23.7 22.34 27.67 20.83 27.57 28.01 24.79 28.96 25.77 28.57 23.52 21.86 28.95 24.53 27.65 27.29 29.07 32.47 29.65 22.11 22.43 20.04]';
>> y=[15.02 14.86 13.98 15.91 15.80 14.32 13.76 15.18 14.20 15.40 15.94 14.33 15.11 13.81 15.58 15.85 15.28 16.40 15.02 15.73 14.75 14.35]';
>> x=[ones(22,1),x1];
>> [b,bint,r,rint,stats]=regress(y,x,0.05)
b =
13.5266
0.0581
bint =
10.8647 16.1886
-0.0442 0.1604
r =
0.1161
0.0352
-1.1546
1.1729
0.6713
-0.8343
-1.2072
-0.0295
-0.8241
0.2131
1.0466
-0.4669
-0.0989
-1.1421
0.4466
0.7375
0.0641
0.9865
-0.2296
0.9185
-0.0801
-0.3412
rint =
-1.4473 1.6796
-1.5023 1.5726
-2.6219 0.3128
-0.2085 2.5544
-0.8647 2.2072
-2.3455 0.6769
-2.6740 0.2596
-1.5745 1.5154
-2.3542 0.7059
-1.3367 1.7630
-0.4326 2.5258
-1.9752 1.0413
-1.6435 1.4456
-2.6185 0.3343
-1.1069 2.0002
-0.7947 2.2698
-1.4781 1.6063
-0.3494 2.3224
-1.7534 1.2941
-0.5486 2.3857
-1.6193 1.4592
-1.7914 1.1091
stats =
0.0656 1.4035 0.2500 0.5710
>> rcoplot(r,rint)
参数回归结果为对应的置信区间分别为[-1.2103 3.2561];;[ -0.0318 0.0824];[1.5526 2.4360];[1.0444 2.7813]
r2=0.9226(越接近于1，回归效果越显著)，F= 63.5671， p=0.0000，由p<0.05, 可知
回归模型
y= 1.0229-0.0253*x1+1.9943*x2+1.9128*x3成立。

再分别求y关于单个变量x1，x2, x3的线性回归方程
y=13.5989+0.0547*x1
y= 5.6018+1.8453*x2
y= 7.1238+4.7430*x3
分别求y关于两个变量x1，x2, x3的线性回归方程
y=1.0113+0.0252*x1+2.6057*x2 y=2.0661+0.0318*x2+2.3961*x3 y=6.0781+0.0806*x1+4.0719*x3
三、优化理论中的线性规划问题---生产安排
1、问题提出：某公司打算利用具有下列成分（见下表）的合金配制一种新型合金100公斤，新合金含铅，锌，锡的比例为3：2：5。

合金品种 1 2 3 4 5 含铅% 含锌% 含锡% 30 60 10
10 20 70
50 20 30
10 10 80
50 10 40 单价（元/k g ）
8.6 6.0 8.9 5.7
8.8
要求：（1）根据题意，列出该问题的线性规划模型；
（2）利用单纯形法求解（1）中的模型，并写出分配方案； （3）编程实现上述求解过程；
（4）利用程序验证上述模型的最优解。

2、程序设计
设每种合金品种取值i x 千克 （12345i =、
、
、、） 根据题意建立线性规划方程得：
目标费用最小
()f x =8.6*x16*x28.9*x3 5.7*x48.8*x5;++++
x1x2x3x4x5100;++++=
(30*x1+10*x2+50*x3+10*x4+50*x5)(60*x1+20*x2+20*x3+10*x4+10*x5)=1.5;
÷()()30*x110*x250*x310*x450*x510*x170*x230*x380*x440*x50.6
++++÷++++=()()60*x120*x220*x310*x410*x510*x170*x230*x380*x440*x50.4
++++÷++++=
结果为：
结果分析
当x1=11.1111 x2=0 x3=44.44444 x4=44.44444 x5=0时
取得费用最小值为 744.4444元
当铅减少1单位时总费用将减少7.4444元 当锌减少1单位时总费用将增加1.703704元 当锡减少1单位时总费用将减少230.0926元
四、非线性方程求解
1、问题提出：分别用二分法、牛顿切线法、迭代法求解非线性方程sin 03
x
x -=的非负实数根。

要求：（1）精确到6
10-，取不同的初值计算，输出初值、根的近似值和迭代次数，分
析根的收敛域。

（2）编写二分法、牛顿切线法的程序。

（可以用Matlab 或C 语言）。

（3）迭代法求解（可构造不同的迭代公式，如12tan n n x x +=等）。

（4）比较三种方法的优劣。

2、模型程序设计求解：
迭代法求解
首先要确定方程实数根存在的大致范围。

为此，先将方程变成标准形式f(x)
=sin 03
x
x -=。

作f(x)的曲线图：
>> x=-2*pi:0.1:2*pi; >> f=sin(x)-x./3; >> plot(x,f);grid on;
从曲线上可以看出，函数的零点大约在x 1 ≈0和x 2 ≈ 2.2附近
(1)直接使用指令fzero 求出方程在x1≈ 0时的根。

>> x1=fzero('sin(x)-x/3',0)
x1 =
(2)若键入：fzero (' in(x)-x/3',0, optimset('disp', 'iter'))，将显示迭代过程。

matlab运行结果为：
中间数据表明，求根过程中不断缩小探测范围，最后得出0满足精度的近似根。

(3)求x2≈ 2.2的根：
(4)
>> x2=fzero (' sin(x)-x/3',2.2, optimset('disp', 'iter'))
结果为：
Search for an interval around 2.2 containing a sign change:
Func-count a f(a) b f(b) Procedure
1 2.
2 0.0751631 2.2 0.0751631 initial
interval
3 2.13777 0.130936 2.26223 0.016260
4 search
5 2.112 0.153089 2.288 -0.00902002 search
Search for a zero in the interval [2.112, 2.288]:
Func-count x f(x) Procedure
5 2.288 -0.00902002 initial
6 2.27821 0.00064474 interpolation
7 2.27886 2.26228e-006 interpolation
8 2.27886 -5.49827e-012 interpolation
9 2.27886 -1.11022e-016 interpolation
10 2.27886 -1.11022e-016 interpolation
Zero found in the interval [2.112, 2.288]
x2 =
2.2789
中间数据表明，求根过程中不断缩小探测范围，最后得出2.2789满足精度的近似根。

二分法求解：
程序为：
建立M文件：
function root=HalfInterval(f,a,b,eps)
if(nargin==3)
eps=1.0e-6;
end
f1=subs(sym(f),findsym(sym(f)),a);
f2=subs(sym(f),findsym(sym(f)),b);
if(f1==0)
root=a;
end
if(f2==0)
root=b;
end
if(f1*f2>0)
disp('两端点函数值乘积大于0! ');
return;
else
root=FindRoots(f,a,b,eps);
end
function r=FindRoots(f,a,b,eps)
f_1=subs(sym(f),findsym(sym(f)),a);
f_2=subs(sym(f),findsym(sym(f)),b);
mf=subs(sym(f),findsym(sym(f)),(a+b)/2);
if (f_1*mf>0) t=(a+b)/2;
r=FindRoots(f,t,b,eps); else
if (f_1*mf==0) r=(a+b)/2; else
if (abs(b-a)<=eps) r=(b+3*a)/4; else
s=(a+b)/2;
r=FindRoots(f,a,s,eps); end
end
在窗口运行： >> f=inline('sin(x)-x/3') HalfInterval(f,1,3) f =
Inline function: f(x) = sin(x)-x/3 ans =
2.2789
x2=2.2789,x0=0;
五、非线性回归问题-------多项式回归
1、问题提出：给动物口服某种药物A 1000mg ，每间隔1小时测定血药浓度（g/ml ），得到表9-5的数据（血药浓度为5头供试动物的平均值）。

血药浓度与服药时间测定结果表：
服药时间x （小时） 1 2 3 4
5 6 7 8 9 血药浓度y (g/ml )
21.89
47.13
61.86
70.78
72.81
66.36
50.34
25.31
3.17
要求：1）画出散点图y 与x ，并观察y 与x 的关系；
2）求y 关于x 的一元线性回归方程：11^
0^
x a a y +=-----（1），求出01,a a 的值； 3）对上述回归模型和回归系数进行检验；
4）再求y 关于x 的一元多项式线性回归方程。

（如： 2
2211^
0^
x a x a a y ++=----（2））求出123,,a a a 的值，并比较二个回归方程对原来问题求解的优劣。

5）编程实现上述求解过程。

2、程序设计：
>> x=[1 2 3 4 5 6 7 8 9];
>> y=[21.89 47.13 61.86 70.78 72.81 66.36 50.34 25.31 3.17]; >> plot(x,y,'b*');title('血药浓度与服药时间'); >> A=[ones(size(x));x]'; >> a=A\y' a =
60.6111 -2.7967
>> b=[ -2.7967 60.6111]; >> y=poly2str(b,'x')
结果为：
由图上可观察出y 与x 之间曲线大致为抛物线。

y =
-2.7967 x + 60.6111
即6111.600=a ，7967.21-=a 。

第二中线性回归方程： >> c=A\y'
c =
-8.3655
34.8269
-3.7624
>> d=[-3.7624 34.8269 -8.3655]; >> y=poly2str(d,'x') 结果为： y =
-3.7624 x^2 + 34.8269 x - 8.3655
即7624.33-=a ，8269.342=a ，3655.80-=a 。

对回归模型和回归系数进行检验：
>> x=[1 2 3 4 5 6 7 8 9]';
>> y=[21.89 47.13 61.86 70.78 72.81 66.36 50.34 25.31 3.17]'; >> x1=[ones(9,1), x];
>> [b,bint,r,rint,stats]=regress(y,x1,0.05) 结果为： b =
60.6111 -2.7967
bint =
17.5913 103.6309 -10.4415 4.8482 r =
-35.9244 -7.8878 9.6389 21.3556 26.1822 22.5289 9.3056 -12.9278 -32.2711
rint =
-72.5693 0.7204
-62.3381 46.5626
-47.6076 66.8854
-34.7095 77.4206
-28.5685 80.9330
-33.1065 78.1643
-47.9923 66.6034
-66.4728 40.6173
-71.9576 7.4154
stats =
0.0966 0.7483 0.4157 627.1365
>> rcoplot(r,rint)
由图中可看出与所求的线性方程非常符合。

数据的残差离零点距离较近，说明回归模型
y =-2.7967 x + 60.6111能较好的符合原始数据。

六、非线性规划问题
1、问题提出：现有两种原料A和B数量分别为1200千克和1500千克,需要分配用于生产3种产品.其中每种产品生产的产量Q与两种原料的关系分别为:
y x y x Q 21005.0),(=,22008.0),(xy y x Q =,xy y x Q 01.0),(3=
每种产品的利润函数为:
2
01.010)(Q Q Q R -=
问:应如何分配,才能使生产三种产品的总利润最大.
要求: 1)介绍非线性规划理论; 2)求出最优解. 2、问题建模： 决策变量：
设 A 数量分给三种产品的数量为321,,x x x ，B 数量分别为321,,y y y ，三种产品的产量分别为1Q ，2Q ，3Q ，则总利润为R 。

目标函数：
12
11005.0y x Q =，2222008.0y x Q =，33301.0y x Q =
2
333322
222222121121)(01.001.010)(01.0008.010)(01.0005.010max y x y x y x y x y x y x R -⨯+-⨯+-⨯=
..t S
1200321≤++x x x 1500321≤++y y y
利用LINGO 软件求得最优解。

结果为：
由上述可知最大的利润为：4725.0 R 。