原子距离问题

实验7 无约束优化
题目5
某分子由25个原子组成，并且已经通过实验测量得到了其中某些原子对之间的距离（假设在平面结构上讨论），如表7.8所示。

请你确定每个原子的位置关系。

【模型建立】
设第i 个点所在的位置为()i i x ,y ，因为所求的是各原子间的位置关系，可以设定第一个点坐标为)0,0(),(11=y x ，然后再计算其他原子的位置()i i x ,y ，使它在最大程度上满足上表中提供的数据，即让
22[()+()-]2
2
i j i j ij i,j
z x -x y -y d =∑
达到最小，其中ij d 表示第i 个原子和第j 个原子之间的距离，数据如上表中所示。

问题转化为无约束优化：
∑--+-j
i ij j i j i d y y x x ,2
222]
)()[(min
【算法设计】
上面是一个无约束优化问题，要求
∑--+-j
i ij j i j i d y y x x ,2
222]
)()[(min
可以归结到无约束规划模型min f(x)，令X=[0,X2,X3,X4,……,x25,0,y2,y3,……,y25],调用基本命
令fminunc来做，也可以令
⎥
⎦
⎤
⎢
⎣
⎡
=
25
24
3
2
25
24
3
2
y
y
y
y
x
x
x
x
x
Λ
Λ
,调用lsqnonlin命令来做。

【程序】
方法一用lsqnonlin命令实现
function f=distance(x,d)
f(1)=(x(1,3))^2+(x(2,3))^2-d(1)^2; %对应条件第4个和第1个原子间距f(2)=(x(1,11))^2+(x(2,11))^2-d(2)^2;
f(3)=(x(1,12))^2+(x(2,12))^2-d(3)^2;
f(4)=(x(1,16))^2+(x(2,16))^2-d(4)^2;
f(5)=(x(1,20))^2+(x(2,20))^2-d(5)^2;
f(6)=(x(1,4)-x(1,1))^2+(x(2，4)-x(2,1))^2-d(6)^2;
f(7)=(x(1,15)-x(1,1))^2+(x(2,15)-x(2,1))^2-d(7)^2;
f(8)=(x(1,16)-x(1,1))^2+(x(2,16)-x(2,1))^2-d(8)^2;
f(9)=(x(1,24)-x(1,1))^2+(x(2,24)-x(2,1))^2-d(9)^2;
f(10)=(x(1,4)-x(1,2))^2+(x(2,4)-x(2,2))^2-d(10)^2;
f(11)=(x(1,19)-x(1,2))^2+(x(2,19)-x(2,2))^2-d(11)^2;
f(12)=(x(1,20)-x(1,2))^2+(x(2,20)-x(2,2))^2-d(12)^2;
f(13)=(x(1,23)-x(1,2))^2+(x(2,23)-x(2,2))^2-d(13)^2;
f(14)=(x(1,4)-x(1,3))^2+(x(2,4)-x(2,3))^2-d(14)^2;
f(15)=(x(1,11)-x(1,3))^2+(x(2,11)-x(2,3))^2-d(15)^2;
f(16)=(x(1,23)-x(1,3))^2+(x(2,23)-x(2,3))^2-d(16)^2;
f(17)=(x(1,7)-x(1,5))^2+(x(2,7)-x(2,5))^2-d(17)^2;
f(18)=(x(1,12)-x(1,5))^2+(x(2,12)-x(2,5))^2-d(18)^2;
f(19)=(x(1,18)-x(1,5))^2+(x(2,18)-x(2,5))^2-d(19)^2;
f(20)=(x(1,24)-x(1,5))^2+(x(2,24)-x(2,5))^2-d(20)^2;
f(21)=(x(1,7)-x(1,6))^2+(x(2,7)-x(2,6))^2-d(21)^2;
f(22)=(x(1,13)-x(1,6))^2+(x(2,13)-x(2,6))^2-d(22)^2;
f(23)=(x(1,15)-x(1,6))^2+(x(2,15)-x(2,6))^2-d(23)^2;
f(24)=(x(1,19)-x(1,6))^2+(x(2,19)-x(2,6))^2-d(24)^2;
f(25)=(x(1,20)-x(1,6))^2+(x(2,20)-x(2,6))^2-d(25)^2;
f(26)=(x(1,13)-x(1,7))^2+(x(2,13)-x(2,7))^2-d(26)^2;
f(27)=(x(1,17)-x(1,7))^2+(x(2,17)-x(2,7))^2-d(27)^2;
f(28)=(x(1,12)-x(1,8))^2+(x(2,12)-x(2,8))^2-d(28)^2;
f(29)=(x(1,14)-x(1,8))^2+(x(2,14)-x(2,8))^2-d(29)^2;
f(30)=(x(1,21)-x(1,8))^2+(x(2,21)-x(2,8))^2-d(30)^2;
f(31)=(x(1,10)-x(1,9))^2+(x(2,10)-x(2,9))^2-d(31)^2;
f(32)=(x(1,12)-x(1,9))^2+(x(2,12)-x(2,9))^2-d(32)^2;
f(33)=(x(1,18)-x(1,9))^2+(x(2,18)-x(2,9))^2-d(33)^2;
f(34)=(x(1,19)-x(1,9))^2+(x(2,19)-x(2,9))^2-d(34)^2;
f(35)=(x(1,21)-x(1,9))^2+(x(2,21)-x(2,9))^2-d(35)^2;
f(36)=(x(1,17)-x(1,10))^2+(x(2,17)-x(2,10))^2-d(36)^2;
f(37)=(x(1,24)-x(1,10))^2+(x(2,24)-x(2,10))^2-d(37)^2;
f(38)=(x(1,14)-x(1,11))^2+(x(2,14)-x(2,11))^2-d(38)^2;
f(39)=(x(1,16)-x(1,11))^2+(x(2,16)-x(2,11))^2-d(39)^2;
f(40)=(x(1,14)-x(1,12))^2+(x(2,14)-x(2,12))^2-d(40)^2;
f(41)=(x(1,18)-x(1,12))^2+(x(2,18)-x(2,12))^2-d(41)^2;
f(42)=(x(1,14)-x(1,13))^2+(x(2,14)-x(2,13))^2-d(42)^2;
f(43)=(x(1,15)-x(1,13))^2+(x(2,15)-x(2,13))^2-d(43)^2;
f(44)=(x(1,19)-x(1,15))^2+(x(2,19)-x(2,15))^2-d(44)^2;
f(45)=(x(1,22)-x(1,15))^2+(x(2,22)-x(2,15))^2-d(45)^2;
f(46)=(x(1,17)-x(1,16))^2+(x(2,17)-x(2,16))^2-d(46)^2;
f(47)=(x(1,18)-x(1,16))^2+(x(2,18)-x(2,16))^2-d(47)^2;
f(48)=(x(1,19)-x(1,18))^2+(x(2,19)-x(2,18))^2-d(48)^2;
f(49)=(x(1,22)-x(1,18))^2+(x(2,22)-x(2,18))^2-d(49)^2;
f(50)=(x(1,23)-x(1,18))^2+(x(2,23)-x(2,18))^2-d(50)^2;
f(51)=(x(1,22)-x(1,20))^2+(x(2,22)-x(2,20))^2-d(51)^2;
f(52)=(x(1,22)-x(1,21))^2+(x(2,22)-x(2,21))^2-d(52)^2;
clear all
x0=[zeros(1,24);ones(1,24)];
d=[0.9607,0.4399,0.8143,1.3765,1.2722,0.5294,0.6144,0.3766,0.6893,...
0.9488,0.8000,1.1090,1.1432,0.4758,1.3402,0.7006,0.4945,1.0559,...
0.6810,0.3587,0.3351,0.2878,1.1346,0.3870,0.7511,0.4439,0.8363,...
0.3208,0.1574,1.2736,0.5781,0.9254,0.6401,0.2467,0.4727,1.3840,...
0.4366,1.0307,1.3904,0.5725,0.7660,0.4394,1.0952,1.0422,1.8255,...
1.4325,1.0851,0.4995,1.2277,1.1271,0.7060,0.8052]';%设定初值
[x,norms]=lsqnonlin(@distance,x0,[],[],[],d);
p=x'; %p第一行即为第二个原子的横、纵坐标；p第二行为第三个原子的横、纵坐标……
a=[0,x(1,:)];
b=[0,x(2,:)];
plot(a,b,'*'); %画散点图表示出原子的位置
【输出结果】
每个原子的位置如下（即第二个原子的位置为（1.2378,0.0746），第三个原子为（1,6142,1,。

1211），……）：
p =
1.2378 0.0746
1.6142 1.1211
0.5485 0.7778
0.9121 0.4836
0.8546 1.1533
0.7009 0.4533
0.9871 0.6573
0.4694 0.2363
0.8326 1.0437 0.4356 0.6211 -0.0650 -0.4184 0.8165 0.1140 0.7351 0.2500 0.3971 0.4986 1.8229 0.2932 1.3241 -0.3717 1.7657 0.9905 1.3514 0.7020 0.8921 0.7787 0.5040 1.1701 1.3106 1.1878 0.8031 1.8060 0.5290 1.4757 0.8946 0.7093
为了形象直观地表示出各点的位置，画出如下的散点图：
-0.5
00.51 1.52
-0.50
0.5
1
1.5
2
方法二 用fminunc 实现
function f=distance1(x,d)
f=(x(1,3))^2+(x(2,3))^2-d(1)^2;
f=f+(x(1,11))^2+(x(2,11))^2-d(2)^2; f=f+(x(1,12))^2+(x(2,12))^2-d(3)^2;
f=f+(x(1,16))^2+(x(2,16))^2-d(4)^2;
f=f+(x(1,20))^2+(x(2,20))^2-d(5)^2;
f=f+(x(1,4)-x(1,1))^2+(x(2*5-2)-x(2,1))^2-d(6)^2;
f=f+(x(1,15)-x(1,1))^2+(x(2,15)-x(2,1))^2-d(7)^2;
f=f+(x(1,16)-x(1,1))^2+(x(2,16)-x(2,1))^2-d(8)^2;
f=f+(x(1,24)-x(1,1))^2+(x(2,24)-x(2,1))^2-d(9)^2;
f=f+(x(1,4)-x(1,2))^2+(x(2,4)-x(2,2))^2-d(10)^2;
f=f+(x(1,19)-x(1,2))^2+(x(2,19)-x(2,2))^2-d(11)^2; f=f+(x(1,20)-x(1,2))^2+(x(2,20)-x(2,2))^2-d(12)^2; f=f+(x(1,23)-x(1,2))^2+(x(2,23)-x(2,2))^2-d(13)^2; f=f+(x(1,4)-x(1,3))^2+(x(2,4)-x(2,3))^2-d(14)^2;
f=f+(x(1,11)-x(1,3))^2+(x(2,11)-x(2,3))^2-d(15)^2; f=f+(x(1,23)-x(1,3))^2+(x(2,23)-x(2,3))^2-d(16)^2; f=f+(x(1,7)-x(1,5))^2+(x(2,7)-x(2,5))^2-d(17)^2;
f=f+(x(1,12)-x(1,5))^2+(x(2,12)-x(2,5))^2-d(18)^2; f=f+(x(1,18)-x(1,5))^2+(x(2,18)-x(2,5))^2-d(19)^2; f=f+(x(1,24)-x(1,5))^2+(x(2,24)-x(2,5))^2-d(20)^2; f=f+(x(1,7)-x(1,6))^2+(x(2,7)-x(2,6))^2-d(21)^2;
f=f+(x(1,13)-x(1,6))^2+(x(2,13)-x(2,6))^2-d(22)^2; f=f+(x(1,15)-x(1,6))^2+(x(2,15)-x(2,6))^2-d(23)^2; f=f+(x(1,19)-x(1,6))^2+(x(2,19)-x(2,6))^2-d(24)^2; f=f+(x(1,20)-x(1,6))^2+(x(2,20)-x(2,6))^2-d(25)^2; f=f+(x(1,13)-x(1,7))^2+(x(2,13)-x(2,7))^2-d(26)^2; f=f+(x(1,17)-x(1,7))^2+(x(2,17)-x(2,7))^2-d(27)^2; f=f+(x(1,12)-x(1,8))^2+(x(2,12)-x(2,8))^2-d(28)^2; f=f+(x(1,14)-x(1,8))^2+(x(2,14)-x(2,8))^2-d(29)^2; f=f+(x(1,21)-x(1,8))^2+(x(2,21)-x(2,8))^2-d(30)^2; f=f+(x(1,10)-x(1,9))^2+(x(2,10)-x(2,9))^2-d(31)^2; f=f+(x(1,12)-x(1,9))^2+(x(2,12)-x(2,9))^2-d(32)^2; f=f+(x(1,18)-x(1,9))^2+(x(2,18)-x(2,9))^2-d(33)^2; f=f+(x(1,19)-x(1,9))^2+(x(2,19)-x(2,9))^2-d(34)^2; f=f+(x(1,21)-x(1,9))^2+(x(2,21)-x(2,9))^2-d(35)^2; f=f+(x(1,17)-x(1,10))^2+(x(2,17)-x(2,10))^2-d(36)^2; f=f+(x(1,24)-x(1,10))^2+(x(2,24)-x(2,10))^2-d(37)^2; f=f+(x(1,14)-x(1,11))^2+(x(2,14)-x(2,11))^2-d(38)^2; f=f+(x(1,16)-x(1,11))^2+(x(2,16)-x(2,11))^2-d(39)^2; f=f+(x(1,14)-x(1,12))^2+(x(2,14)-x(2,12))^2-d(40)^2; f=f+(x(1,18)-x(1,12))^2+(x(2,18)-x(2,12))^2-d(41)^2; f=f+(x(1,14)-x(1,13))^2+(x(2,14)-x(2,13))^2-d(42)^2; f=f+(x(1,15)-x(1,13))^2+(x(2,15)-x(2,13))^2-d(43)^2; f=f+(x(1,19)-x(1,15))^2+(x(2,19)-x(2,15))^2-d(44)^2; f=f+(x(1,22)-x(1,15))^2+(x(2,22)-x(2,15))^2-d(45)^2; f=f+(x(1,17)-x(1,16))^2+(x(2,17)-x(2,16))^2-d(46)^2; f=f+(x(1,18)-x(1,16))^2+(x(2,18)-x(2,16))^2-d(47)^2;
f=f+(x(1,19)-x(1,18))^2+(x(2,19)-x(2,18))^2-d(48)^2;
f=f+(x(1,22)-x(1,18))^2+(x(2,22)-x(2,18))^2-d(49)^2;
f=f+(x(1,23)-x(1,18))^2+(x(2,23)-x(2,18))^2-d(50)^2;
f=f+(x(1,22)-x(1,20))^2+(x(2,22)-x(2,20))^2-d(51)^2;
f=f+(x(1,22)-x(1,21))^2+(x(2,22)-x(2,21))^2-d(52)^2;
clear all;
d=[0.9607,0.4399,0.8143,1.3765,1.2722,0.5294,0.6144,0.3766,0.6893,0.9488,0.8000,1.1090,1.1432,0 .4758,1.3402,0.7006,0.4945,1.0559,0.6810,0.3587,0.3351,0.2878,1.1346,0.3870,0.7511,0.4439,0.836 3,0.3208,0.1574,1.2736,0.5781,0.9254,0.6401,0.2467,0.4727,1.3840,0.4366,1.0307,1.3904,0.5725,0.7 660,0.4394,1.0952,1.0422,1.8255,1.4325,1.0851,0.4995,1.2277,1.1271,0.7060,0.8052]';
x0=[zeros(1,26),ones(1,24)];
opt = optimset('tolx',1e-16,'tolf',1e-16);
[x,z,exit1,out1] = fminunc(@distance1,x0,opt,d);
p=[0,0;x(1:24)',x(25:48)'];
a=[0,x(1:24)];
b=[0,x(25:48)];
plot(a,b,'*');%画散点图表示原子位置
【输出结果】
每个原子的位置如下（即第二个原子的位置为（-0.6818，0.6856），第二个原子为（-0.3563，0.2298），……）：
p =
0 0
-0.6818 0.6856
-0.3563 0.2298
-0.0530 0.9715
-0.4866 1.1721
0.0506 1.0700
0.3992 0.6744
0.0547 0.5349
-0.3948 0.3789
-0.1388 1.0752
0.4190 0.9660
-0.3796 -0.3227
-0.7430 0.3687
0.1511 0.9167
-0.2615 0.6897
-0.5621 0.0830
-0.9836 0.9375
-0.2070 -0.2659
-0.0144 0.4168
0.0717 0.9085
-0.0117 1.2804
-0.0890 1.6054
0.6877 1.4172
-0.6661 1.3305
0.0125 0.7379
画出位置散点图如下：
-1-0.8-0.6-0.4-0.200.20.40.60.8
【结果分析】
（1）关于lsqononlin与fminunc精确性的探讨
用上面两种不同的命令所求出来的结果并不一样，可能是因为模型建立稍有差别导致算法不同而引起的。

为了找到更适合解决本题的模型，在命令窗口中输入以下命令比较两种算法的精确性。

对于方法一（用lsqnonlin命令实现）
f(2)=(x(1,11))^2+(x(2,11))^2-d(2)^2;
f(3)=(x(1,12))^2+(x(2,12))^2-d(3)^2;
f(4)=(x(1,16))^2+(x(2,16))^2-d(4)^2;
f(5)=(x(1,20))^2+(x(2,20))^2-d(5)^2;
f(6)=(x(1,4)-x(1,1))^2+(x(2*5-2)-x(2,1))^2-d(6)^2;
f(7)=(x(1,15)-x(1,1))^2+(x(2,15)-x(2,1))^2-d(7)^2;
f(8)=(x(1,16)-x(1,1))^2+(x(2,16)-x(2,1))^2-d(8)^2;
f(9)=(x(1,24)-x(1,1))^2+(x(2,24)-x(2,1))^2-d(9)^2;
f(10)=(x(1,4)-x(1,2))^2+(x(2,4)-x(2,2))^2-d(10)^2;
f(11)=(x(1,19)-x(1,2))^2+(x(2,19)-x(2,2))^2-d(11)^2; f(12)=(x(1,20)-x(1,2))^2+(x(2,20)-x(2,2))^2-d(12)^2; f(13)=(x(1,23)-x(1,2))^2+(x(2,23)-x(2,2))^2-d(13)^2; f(14)=(x(1,4)-x(1,3))^2+(x(2,4)-x(2,3))^2-d(14)^2;
f(15)=(x(1,11)-x(1,3))^2+(x(2,11)-x(2,3))^2-d(15)^2; f(16)=(x(1,23)-x(1,3))^2+(x(2,23)-x(2,3))^2-d(16)^2; f(17)=(x(1,7)-x(1,5))^2+(x(2,7)-x(2,5))^2-d(17)^2;
f(18)=(x(1,12)-x(1,5))^2+(x(2,12)-x(2,5))^2-d(18)^2; f(19)=(x(1,18)-x(1,5))^2+(x(2,18)-x(2,5))^2-d(19)^2; f(20)=(x(1,24)-x(1,5))^2+(x(2,24)-x(2,5))^2-d(20)^2; f(21)=(x(1,7)-x(1,6))^2+(x(2,7)-x(2,6))^2-d(21)^2;
f(22)=(x(1,13)-x(1,6))^2+(x(2,13)-x(2,6))^2-d(22)^2; f(23)=(x(1,15)-x(1,6))^2+(x(2,15)-x(2,6))^2-d(23)^2; f(24)=(x(1,19)-x(1,6))^2+(x(2,19)-x(2,6))^2-d(24)^2; f(25)=(x(1,20)-x(1,6))^2+(x(2,20)-x(2,6))^2-d(25)^2; f(26)=(x(1,13)-x(1,7))^2+(x(2,13)-x(2,7))^2-d(26)^2; f(27)=(x(1,17)-x(1,7))^2+(x(2,17)-x(2,7))^2-d(27)^2; f(28)=(x(1,12)-x(1,8))^2+(x(2,12)-x(2,8))^2-d(28)^2; f(29)=(x(1,14)-x(1,8))^2+(x(2,14)-x(2,8))^2-d(29)^2; f(30)=(x(1,21)-x(1,8))^2+(x(2,21)-x(2,8))^2-d(30)^2; f(31)=(x(1,10)-x(1,9))^2+(x(2,10)-x(2,9))^2-d(31)^2; f(32)=(x(1,12)-x(1,9))^2+(x(2,12)-x(2,9))^2-d(32)^2; f(33)=(x(1,18)-x(1,9))^2+(x(2,18)-x(2,9))^2-d(33)^2; f(34)=(x(1,19)-x(1,9))^2+(x(2,19)-x(2,9))^2-d(34)^2; f(35)=(x(1,21)-x(1,9))^2+(x(2,21)-x(2,9))^2-d(35)^2; f(36)=(x(1,17)-x(1,10))^2+(x(2,17)-x(2,10))^2-d(36)^2; f(37)=(x(1,24)-x(1,10))^2+(x(2,24)-x(2,10))^2-d(37)^2; f(38)=(x(1,14)-x(1,11))^2+(x(2,14)-x(2,11))^2-d(38)^2; f(39)=(x(1,16)-x(1,11))^2+(x(2,16)-x(2,11))^2-d(39)^2; f(40)=(x(1,14)-x(1,12))^2+(x(2,14)-x(2,12))^2-d(40)^2; f(41)=(x(1,18)-x(1,12))^2+(x(2,18)-x(2,12))^2-d(41)^2; f(42)=(x(1,14)-x(1,13))^2+(x(2,14)-x(2,13))^2-d(42)^2; f(43)=(x(1,15)-x(1,13))^2+(x(2,15)-x(2,13))^2-d(43)^2; f(44)=(x(1,19)-x(1,15))^2+(x(2,19)-x(2,15))^2-d(44)^2; f(45)=(x(1,22)-x(1,15))^2+(x(2,22)-x(2,15))^2-d(45)^2;
运行结果为
a =0.2190 对于方法二（用fminunc命令实现）
f =0.0375
两种算法初值的设定是一样的。

对比两个运行结果，用fminunc求出的距离差的和小于用lsqonlin求出的距离差的和。

由此可见，用fminunc命令求解更为精确。

（2）改变初值对结果的影响
鉴于前面的分析，我们采用fminunc命令得到的结果来继续下面的讨论。

将初值改为
x0=[zeros(1,2),ones(1,46)];
运行结果如下：
p =
0 0
0.5553 1.5928
0.1561 0.4287
0.7454 0.6616
0.8003 1.1268
1.1061 0.7317
1.1218 0.9560
0.8743 1.1478
0.2187 0.4924
0.9728 0.8388
0.3939 0.8547
-0.4765 0.1503
0.0492 0.8143
0.9413 0.6822
0.5017 0.4869
0.0320 1.2791
0.2655 1.3360
1.6969 1.3201
0.5989 0.2906
0.8971 0.6896
1.2488 0.2158
1.2155 1.2784
1.7608 0.6915
0.7832 1.3882
0.8195 0.9554画出位置散点图：
对于其他初值的设定情况也大致相同。

由此可见，改变初值，对于最终结果的影响是很大的。

不同的初值，得到原子的位置不一样。

我觉得可能的原因是fminunc 求出的是局部极小值，在不同的初值条件下函数会收敛到不同的极小值，因此会有上述情况发生。

【实验结论】
对于无约束优化求极小值的问题，可以采用lsqnonlin 和fminunc 等不同的命令，得到的结果虽算法的不同而有所不同。

此外，设置不同的初值，也可以得到不同的结果，因此说明初值的选取对研究问题其一定的作用。

为了得到更准确的模型，需要根据实际问题慎重地选取初值。

题目8 给药方案设计需要依据药
物吸收与排除过程的原理。

药物进入机体后随血液输送到全身，不断地被吸收、分布、代谢，最终排出体外。

药物在血液中的浓度，即单
位体积血液中的药物含量，称血药浓度。

在最简单的一室模型中，将
整个机体看作一个房室，称中心室，室内的血药浓度是均匀的。

这里我们用一室模型，讨论在口服给药方式下血药浓度的变化规律，及根据实验数据拟合参数的方法。

吸收室c1（t ）
中心室 c （t ）
k1
k
口服给药方式相当于先有一个将药从肠胃吸收入血液的过程，这个过程可简化为在药物进入中心室之前有一个吸收室（如图），记中心室和吸收室的容积分别为V ，V1，而t 时刻的血药浓度分别为c （t ），c1（t ）；中心室的排除速率为k ，吸收速率为k1（这里k 和k1分别是中心室和吸收室血药浓度变化率与浓度本身的比例系数）。

设t=0时刻口服剂量为d 的药物，容易写出吸收室的血药浓度c1（t ）的微分方程为
1
1111,(0).dc k c dt d c V ⎧=-⎪⎪⎨⎪=⎪⎩
中心室血药浓度c （t ）的变化率由两部分组成：与c 成正比的排除（比例系数ｋ）；与ｃ１成正比的吸收（比例系数ｋ１）。

再考虑到中心室和吸收室的容积分别为Ｖ，Ｖ１，得到ｃ（ｔ）的微分方程为
111,(0)0.
V dc kc k c dt V c ⎧=-+⎪
⎨⎪=⎩
由以上两个微分方程不难解出中心室血药浓度
111()().k t kt k d c t e e V k k
--=--
在制定给药方案时必须知道这种药物的3个参数k ，k1，b(=d/v)，实际中通常通过实验数据确定。

设t=0时刻口服一定剂量的药物，下表是实验数据c(t)，请由此确定k ，k1，b 。

【模型建立】
问题可以转换为利用数据拟合方程
)()(111
t k kt e e k
k k b
t c ----=
b=d/V
这个方程是非线性的。

由于这是一个非线性最小二乘拟合的问题，故可以选用MATLAB中的Isqnonlin命令或Isqcurvefit命令。

【程序】
function f=medicine(x,t)%x(1)--b,x(2)--k1,x(3)--k
f = x(1)*x(2)/(x(2) - x(3))*(exp(-x(3)*t)-exp(-x(2)*t));
方法一用Isqcurvefit命令
x0=[1,1,0];
t = [0.083 0.167 0.25 0.50 0.75 1.0 1.5 2.25 3.0 4.0 6.0 8.0 10.0 12.0];
c = [10.9 21.1 27.3 36.4 35.5 38.4 34.8 24.2 23.6 15.7 8.2 8.3 2.2 1.8];
opt = optimset;
[x,norm,res,ef,out,lam] = lsqcurvefit(@medicine,x0,t,c,[],[],opt);
【输出结果】
方法二用Isqnonlin命令
x0 = [50,1,0];
t = [0.083 0.167 0.25 0.50 0.75 1.0 1.5 2.25 3.0 4.0 6.0 8.0 10.0 12.0];
c = [10.9 21.1 27.3 36.4 35.5 38.4 34.8 24.2 23.6 15.7 8.2 8.3 2.2 1.8];
opt=optimset('LargeScale','on','MaxFunEvals',1000,'MaxIter',200);
opt=optimset(opt,'tolx',1e-16,'tolf',1e-16);
[x,norm,res,ef,out,lam]=lsqnonlin(@medicine,x0,[],[],[],t,c);
【输出结果】
本上是一致的。

【结果分析】
（1）Isqnonlin和Isqcurvefit的对比
可以分别使用Isqnonlin和Isqcurvefit求解非线性最小二乘拟合。

实验结果发现结果基本一致。

这是由两者所用的算法完全一样所决定的。

此外，还可以对比分析采用分析导数（雅克比矩阵）、对比分析步长一维搜索算法（混合二三次插值和三次插值）等对Matlab的非线性二乘拟合算法进行探索，所得结果都是一致的。

（2）LM法和GN法的对比
x0=[1,1,0];
t = [0.083 0.167 0.25 0.50 0.75 1.0 1.5 2.25 3.0 4.0 6.0 8.0 10.0 12.0];
c = [10.9 21.1 27.3 36.4 35.5 38.4 34.8 24.2 23.6 15.7 8.2 8.3 2.2 1.8];
opt1=optimset('MaxFunEvals',2000);
[x1,norm1,res1]=lsqcurvefit(@medicine,x0,t,c,[],[],opt1); % LM法
opt2=optimset(opt1,'LevenbergMarquardt','off');
[x2,norm2,res2]=lsqcurvefit(@medicine,x0,t,c,[],[],opt2); % GN法
运行结果如下：
lsqcurvefit函数调用时默认使用的下降方向算法为LM法。

可以通过更改控制参数中的LevenbergMarquard而使用GN算法求解非线性最小二乘法。

结果发现，两种下降方向算法所得结果是完全一致的。

【实验结论】
利用lsqnonlin命令和lsqcurvefit命令作数据拟合，得到的结果基本一致。

对于本题，可以得到b =46.8275 ,k1=3.6212 ,k=0.2803。

用Isqnonlin命令所得到的结果要比用Isqcurvefit命令所得到的结果稍微精确一点。

但是用Isqnonlin命令对初值的设定比较敏感，而用Isqcurvefit命令在较大的初值范围内都能保持比较精确的计算。

此外，利用本题模型对比分析LM法和GN法，分析采用分析导数（雅克比矩阵）、对比分析步长一维搜索算法（混合二三次插值和三次插值）等。

对于本题的模型而言，所得结果都是一致的。