数学实验8月13日作业

1.取不同的初值计算下列平方和形式的非线性规划，尽可能求出所有局部极小点，进 而找出全局极小点，并对不同算法（搜索方向、搜索步长、数值梯度与分析梯度等）的结 果进行分析、比较。

(2). (
)(
)
2
2
2
22
121212min 12114949812324681x x x x x x +-++++-,
(4).()()212222
23
12123min10010,1x x x x x x θ⎧⎫⎡⎤⎪⎪-++-+⎡⎤⎨⎬⎢⎥⎣⎦⎣⎦⎪⎪⎩⎭
，其中 ()()()21112211
1
arc ,02,11arc ,0
22tg x x x x x tg x x x π
θπ
⎧>⎪⎪=⎨⎪+<⎪⎩.
解：(2)由题意设()()(
)
22
2
22
121
21
212,12114949812324681f x x x x x x x x =+-++++-
先求12,f f x x δδδδ（文件名为a.m ），再求1
2
00f x f x δδδδ⎧=⎪⎪⎨⎪=⎪⎩（文件名为b.m ）,解得近似解为12 2.900.22x x =-⎧⎨
=⎩或12 3.460.08x x =⎧⎨
=-⎩或1220.6234.53x x =⎧⎨=-⎩或1221.1836.47x x =-⎧⎨=-⎩或120.280.28x x =⎧⎨=⎩或1222.37
17.49
x x =⎧⎨=-⎩或
1224.1317.92x x =-⎧⎨
=-⎩或120.8323.72x x =-⎧⎨=-⎩或12
0.17
47.72x x =-⎧⎨=-⎩,取初值()0 2.90,0.22x =-（exam0705grad_run.m ）利用matlab 计算可得以下表中数据：
当取初值()0 3.46,0.08x =-（exam07051grad_run.m ）利用matlab 计算可得以下表中数据：
取不同初值求解（只需在程序中不断改变初值即可）。

当初值取()021.18,36.47x =--时可得到以下数据
综上全局极小点为
()
-21.181,-36.468，局部极小点
()
-2.8996,0.21584，
()-21.181,-36.468。

（2）
8． 取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法（搜索方向、搜索步长、数值梯度与分析梯度等）的结果进行分析、比较。

(2)111222
11
min ()()()()T T
z x a x a c x a x a c =-
---+--+,2x ∈ℜ 其中c = (0.7,0.73)，T 1a = (4,4) ,T 2 a = (2.5,3.8)。

解：(2)
1. ①设计程序如exam070702fun.m.建立功能函数，便于调用。

function [f,g]=exam070702fun(x) a1=[4,4]'; a2=[
2.5,
3.8]'; c=[0.7,0.73];
f=1/((x-a1)'*(x-a1)+c(1))-1/((x-a2)'*(x-a2)+c(2)); if nargout>1
g(1)=2*(x(1)-a1(1))/((x(1)-a1(1))^2+(x(2)-a1(2))^2+c(1))^2+2*(x(1)-a2(1))/((x(1)-a2(1))^2+(x(2)-a2(2))^2+c(2))^2 %该为梯度函数的两个方向 g(2)=2*(x(2)-a1(2))/((x(1)-a1(1))^2+(x(2)-a1(2))^2+c(1))^2+2*(x(2)-a2(2))/((x(1)-a2(1))^2+(x(2)-a2(2))^2+c(2))^2 end
②设计脚本程序如exam070702.m ，
opt1=optimset('LargeScale','off','MaxfunEvals',1000,'Tolfun',1e-8,'To lX',1e-8,'Gradobj','on');
opt2=optimset(opt1,'HessUpdate','dfp');
opt3=optimset(opt1,'HessUpdate','steepdesc');
opt4=optimset(opt1,'lineSearchType','cubicpoly'); opt5=optimset(opt5,'HessUpdate','dfp');
opt6=optimset(opt5,'HessUpdate','steepdesc');
%x0=[3,3]';
%x0=[4,4]';
%x0=[5,5]';
%x0=[2,7]';
%x0=[10,10]';
%x0=[50,20]';
%x0=[1,10]';
%x0=[50,100]';
%x0=[43,2]';
x0=[20,43]';
[x1,v1,exit1,out1]=fminunc(@exam070702fun,x0,opt1) [x2,v2,exit2,out2]=fminunc(@exam070702fun,x0,opt2)
[x3,v3,exit3,out3]=fminunc(@exam070702fun,x0,opt3)
[x4,v4,exit4,out4]=fminunc(@exam070702fun,x0,opt4)
[x5,v5,exit5,out5]=fminunc(@exam070702fun,x0,opt5)
[x6,v6,exit6,out6]=fminunc(@exam070702fun,x0,opt6)
③运行结果及分析。

得出答案；制成表格
初值设定为（20，43）
可以看到采用分析方法(给定梯度函数) 比起默认的数值方法得到结果所需要的调用次数明显减少。

DFP算法的函数调用次数较少。

改变初值，采用默认搜索方向BFGS，默认搜索步长，给定梯度函数，得到的结果
结论：由上面计算结果得出，无论采用那种计算方法最终都可以得到函数的极小值的，而采用DFP算法函数的调用次数较少；分析法计算梯度函数效果较好；最终结果取哪一个局部极值同初值密切相关。

5.某分子由25 个原子组成，并且已经通过实验测量得到了其中某些原子对之间的距离
解:每个原子的位置都是未知的，在坐标系中只有相对的位置参数，不妨固定原子1的坐标为(0,0)。

并且分子可以在平面内任意旋转，很难确定每个原子的绝对位置，表中给出了52组数据，可以得到52个方程，而未知数个数为48个（每个原子的x ，y 坐标）。

属于超定方程组，没有确定的解，只能求得最优解。

故采用最小方差标准来求最小二乘解得到最优解。

②模型建立
设第i 个点所在的位置为(,)i i x y ，因为所求的是各原子间的位置关系，可以设定第一个点坐标为11(,)(0,0)x y =，然后再计算其他原子的位置(,)i i x y ，使它在最大
程度上满足上表中提供的数据，即让2
222
10,z=()()i i ij i j
x x y y d ⎡⎤-+--⎣⎦∑达到最小，
求得的坐标值就是最优的解。

其中ij d 表示第i 个原子和第j 个原子之间的距离，
数据如上表所示，问题转化为无约束优化；2
222
10,min=()()i i ij i j
x x y y d ⎡⎤-+--⎣⎦∑
③算法设计
根据题目中所给的数据，利用matlab 编写函数文件（文件名为diatance.m ）,再 主程序为exam5.m, 输出结果
每个原子的位置如下（即第二个原子的位置为（1.2378,0.0746），第三个原子为（1,6142,1,。

1211），……） p =
0.3339 -1.2619 1.5092 -1.3714 0.8438 -0.4442 0.6956 -0.8857
1.3388 -0.6804 0.6158 -0.6627 0.9014 -0.9378 0.3552 -0.4180 1.2183 -0.6526 0.6709 -0.4580 -0.4305 -0.0207 0.2989 -0.7737 0.5112 -0.7080 0.5654 -0.2765 0.6723 -1.7884 -0.1065 -1.3723 1.2895 -1.6922 0.9565 -1.2228 0.9735 -0.7787 1.2460 -0.2878 1.4318 -1.0946 1.9280 -0.4652 1.5063 -0.2325 0.9215 -0.8439
为了形象直观地表示出各点的位置，画出如下的散点图
-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-0.207. 经济学中著名的柯布—道格拉斯（Cobb —Douglas ）生产函数的一般形式为
(,)Q K L aK L αβ=，0α<，1β>
其中Q ,K ,L 分别表示产值、资金、劳动力，式中α,β,a 要由经济统计数据确定。

现有《中国统计年鉴（2003）》给出的统计数据如下表，请用非线性最小二乘拟合求出式中的α，β，a ，并解释α，β的含义。

业人员”。

解：①模型分析
题目已经确定了此数据的拟合形式为Cobb-Douglas 生产函数方程，再以这些数据分别给出线性与非线性的最小二乘拟合，则确定式中的a ,β,α。

在构造函数时，以x 向量的分量分别表示a ,α,β,以t 矩阵的第一列表示资金K ，第二列表示劳动力L 。

②非线性拟合程序代码
先用mutlab 编写函数（文件名为examfun7.m ）,主程序为exam7.m,输出结果如下 x =
0.8345 0.7736 0.7312
norm =
2.7489
res =
Columns 1 through 8
-0.3048 0.0404 0.1000 0.1339 0.1254 -0.1437 -0.0866 -0.0577
Columns 9 through 16
0.1528 0.6527 0.4133 -0.0455 -0.2861 -0.4392 -0.0194 -0.0217
Columns 17 through 19
-1.0567 0.1571 0.7350 ef =
3
out =
firstorderopt: 4.7566e-004 iterations: 14 funcCount: 60 cgiterations: 0
algorithm: 'large-scale: trust-region reflective Newton' message: [1x460 char]
非线性拟合给出的结果，a=0.83444,=0.7736α ,=0.73117β，其误差平方和
norm=2.7489，稍大。

在接下来的分析中将看到它与原数据的拟合程度。

③线性拟合程序代码 1. 作图观察QKL 的趋势
用mutlab 编写程序（文件名为exam71.m ）,运行，得到如下图像
1984
198619881990199219941996199820002002
02
4
6
8
10
12
2. 对数处理 从形式上观察，我们不妨将生产函数适当处理，将所给的Cobb-Douglas 方程两边取对数，则可将非线性方程化为线性形式：ln (,)ln ln ln Q K L a K L αβ=++ 再令ln y Q =,ln k K =,ln l L =ln m a =,则方程可以转化为y m k l αβ=++。

令1
121919111k l k
l k l φ⎛⎫ ⎪
⎪
= ⎪
⎪⎝⎭
，(,,)T x a αβ=，12319(,,,)v v v v v =
则问题转化为超定方程x v φ=
根据最小二乘准则计算可以得到矩阵中的待定系数。

由于要求线性拟合，故采用lsqlin 函数。

3. 线性拟合代码
用matlab 编写线性拟合代码（文件名为exam72.m ），输出结果如下，
x = -0.0066453 0.73489 0.64329 resnorm = 0.35783 residual = Columns 1 through 4 0.42453 -0.13017 -0.16924 -0.17176 Columns 5 through 8 -0.13504 0.033093 0.0068534 -0.0045973 Columns 9 through 12
-0.072331 -0.16972 -0.070651 0.028884 Columns 13 through 16 0.070479 0.092607 0.040343 0.043375 Columns 17 through 19 0.15767 0.035332 -0.0096478 exit = 3 out
= iterations: 6 algorithm: [1x43 char] firstorderopt: 0.0114 cgiterations: 5 message: [1x123 char] a = 0.99338 alpha = 0.73489 beta = 0.64329 4.直接左除
程序名为exam73.m,代码如下：
k=log(K);
l=log(L);
f=log(Q);
M=[ ones(19,1),k,l];
s=M\f
a=exp(s(1))
alpha=s(2)
beta=s(3)
输出结果为：s =
-1.1464
0.6605
1.2627
a =
0.3178
alpha =
0.6605
beta =
1.2627
不符合实际情况，应是左除无法加入参数条件限制导致的。

下对此情况不予讨论4.作图对比
程序为exam74.m,代码如下：for t=1:19
Q1(t)=0.83444*K(t)^0.7736*L(t)^0.73117;
end
for i=1:19
Q2(i)=0.99338*K(i)^0.73489*L(i)^0.64329;
end
ts=1984:2002;
plot(ts,Q,'k'), hold on
plot(ts,Q1,'r'), hold on
plot(ts,Q2,'g') title('原始数据与拟合函数对比')
1984198619881990199219941996199820002002
02
4
6
8
10
12
原始数据与拟合函数对比
黑线为原始数据的折线，红色为非线性拟合，绿色为取对数后的线性拟合。

由于原始数据只有19个，量比较小，还不足以看出拟合的效果究竟如何。

从图上看，它们都基本描述出了Q 的走势，理论上讲，线性拟合的效果更好一些。

分析，对数线性拟合的norm=0.3579，比起用非线性拟合的数据norm=2.4789，此norm 值要小得多，这说明用取对数后用线性拟合效果很好。

而且在试探的过程中，发现结果对给定参数的范围的变化较敏感。

但是对初值的变化很不敏感。

6. 参数含义
通过对柯布-道格拉斯生产函数传递变形后，进行求解得出α,β的值，同样也进行预测数据和原始数据比较。

模型中参数α,β的解释：α是劳动力产出的弹性系数，β是资本产出的弹性系数，从这个模型看出，决定工业系统发展水平的主要因素是投入的劳动力数、固定资产和综合技术水平（包括经营管理水平、劳动力素质、引进先进技术等）.根据α和β的组合情况，它有三种类型：①1
αβ+>
称为递增报酬型，表明按现有技术用扩大生产规模来增加产出是有利的.②
1αβ+<称为递减报酬型，
表明按现有技术用扩大生产规模来增加产出是得不偿失的.③1αβ+=称为不变报酬型，表明生产效率并不会随着生产规模的扩大而提高，只有提高技术水平，才会提高经济效益.。