数学实验作业8月10日第五组

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

2)min(x1212x211)2(49x1249x2281x12324x2681)2 4)min100

1arctg(), >0 2其中(,)

11arctg(), <0 22解：2）在matlab中输入程序（exam0705grad.m，exam0705grad_run.m）：

function [f,g]=exam0705grad(x)

f=(x(1)^2+12*x(2)-11)^2+(49*x(1)^2+49*x(2)^2+81*x(1)+2324*x(2)-681)^2;% compute the function value at x

if nargout > 1 % fun called with 2 output arguments

g(1)=2*(98*x(1) + 81)*(49*x(1)^2 + 81*x(1) + 49*x(2)^2 + 2324*x(2) - 681) + 4*x(1)*(x(1)^2 + 12*x(2) - 11); % compute the gradient evaluated at x

g(2)=288*x(2) + 2*(98*x(2) + 2324)*(49*x(1)^2 + 81*x(1) + 49*x(2)^2 + 2324*x(2) - 681) + 24*x(1)^2 - 264; end clear all; clc;

% comparing different algorithms: using gradient vector format short e x0=[-1.9,2];

'---case1: bfgs, hybrid 2,3 poly-------'

opt1=optimset('LargeScale','off', 'MaxFunEvals',1000,'GradObj','on'); [x1,v1,exit1,out1]=fminunc(@exam0705grad,x0,opt1) pause

'---case2: dfp, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','dfp');

[x2,v2,exit2,out2]=fminunc(@exam0705grad,x0,fopt) pause

'---case3: bfgs, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','bfgs');

[x3,v3,exit3,out3]=fminunc(@exam0705grad,x0,fopt) pause

'---case4: steep, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','steepdesc'); [x4,v4,exit4,out4]=fminunc(@exam0705grad,x0,fopt) pause

'---case5: bfgs, 3rd poly-------'

opt2=optimset(opt1,'LineSearchType','cubicpoly'); [x5,v5,exit5,out5]=fminunc(@exam0705grad,x0,opt2) pause

'---case6: dfp, 3rd poly-------' fopt=optimset(opt2,'HessUpdate','dfp');

[x6,v6,exit6,out6]=fminunc(@exam0705grad,x0,fopt) pause

'---case7:bfgs , 3rd poly-------' fopt=optimset(opt2,'HessUpdate','bfgs');

[x7,v7,exit7,out7]=fminunc(@exam0705grad,x0,fopt) pause

'---case8: steep, 3rd poly-------'

fopt=optimset(opt2,'HessUpdate','steepdesc'); [x8,v8,exit8,out8]=fminunc(@exam0705grad,x0,fopt)

pause

'++++ results of solutions ++++++' solutions=[x1;x2;x3;x4;x5;x6;x7;x8]; funvalues=[v1;v2;v3;v4;v5;v6;v7;v8];

iterations=[out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount;out7.funcCount;out8.funcCount]; [solutions,funvalues,iterations]

则输出： x1 =

-2.8997e+000 2.1583e-001 v1 =

3.2329e-006 exit1 = 1 out1 =

iterations: 32 funcCount: 41 stepsize: 1

firstorderopt: 2.3873e+000

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x438 char] x2 =

-2.3709e+000 2.5629e-001 v2 =

6.8300e+000 exit2 = 0 out2 =

iterations: 401 funcCount: 422 stepsize: 1

firstorderopt: 5.7455e+003

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x128 char] x3 =

-2.8997e+000 2.1583e-001

v3 =

3.2329e-006 exit3 = 1 out3 =

iterations: 32 funcCount: 41 stepsize: 1

firstorderopt: 2.3873e+000

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x438 char] x4 =

-2.0594e+000 2.7379e-001 v4 =

1.2063e+001

exit4 = 0 out4 =

iterations: 118 funcCount: 1000 stepsize: 4.7397e-006 firstorderopt: 1.7291e+002

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x144 char] x5 =

-2.8997e+000 2.1583e-001 v5 =

3.2329e-006 exit5 = 1 out5 =

iterations: 32 funcCount: 41 stepsize: 1

firstorderopt: 2.3873e+000

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x438 char] x6 =

-2.3709e+000 2.5629e-001 v6 =

6.8300e+000 exit6 = 0 out6 =

iterations: 401 funcCount: 422 stepsize: 1

firstorderopt: 5.7455e+003

algorithm: 'medium-scale: Quasi-Newton line search'

message: [1x128 char] x7 =

-2.8997e+000 2.1583e-001 v7 =

3.2329e-006 exit7 = 1 out7 =

iterations: 32 funcCount: 41 stepsize: 1

firstorderopt: 2.3873e+000

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x438 char] x8 =

-2.0594e+000 2.7379e-001

v8 =

1.2063e+001 exit8 = 0 out8 =

iterations: 118 funcCount: 1000 stepsize: 4.7397e-006 firstorderopt: 1.7291e+002

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x144 char] ans =

++++ results of solutions ++++++ ans =

-2.8997e+000 2.1583e-001 3.2329e-006 4.1000e+001 -2.3709e+000 2.5629e-001 6.8300e+000 4.2200e+002 -2.8997e+000 2.1583e-001 3.2329e-006 4.1000e+001 -2.0594e+000 2.7379e-001 1.2063e+001 1.0000e+003

-2.8997e+000 2.1583e-001 3.2329e-006 4.1000e+001 -2.3709e+000 2.5629e-001 6.8300e+000 4.2200e+002 -2.8997e+000 2.1583e-001 3.2329e-006 4.1000e+001 -2.0594e+000 2.7379e-001 1.2063e+001 1.0000e+003

改变x0,程序见附序（di1tix0.m），因太过繁琐，只输出最后x的值，输出： ans =

-2.8956e+000 2.1618e-001 4.6001e-004 4.9000e+001 -2.9377e+000 2.1295e-001 1.0877e+000 4.6600e+002 -2.8956e+000 2.1618e-001 4.6001e-004 4.9000e+001 -1.4106e+000 2.9837e-001 2.9482e+001 1.0000e+003 -2.8956e+000 2.1618e-001 4.6001e-004 4.9000e+001 -2.9377e+000 2.1295e-001 1.0877e+000 4.6600e+002 -2.8956e+000 2.1618e-001 4.6001e-004 4.9000e+001 -1.4106e+000 2.9837e-001 2.9482e+001 1.0000e+003 4）程序见（exam07052grad.m， Untitled.m），输出结果：

-1.9000e+000 1.0000e+000 2.0000e+000 5.0213e+002 1.6000e+001 -1.9000e+000 1.0000e+000 2.0000e+000 5.0213e+002 1.6000e+001 -1.9000e+000 1.0000e+000 2.0000e+000 5.0213e+002 1.6000e+001 -1.9000e+000 1.0000e+000 2.0000e+000 5.0213e+002 1.6000e+001 -1.9000e+000 1.0000e+000 2.0000e+000 5.0213e+002 1.6000e+001 -1.9000e+000 1.0000e+000 2.0000e+000 5.0213e+002 1.6000e+001 -1.9000e+000 1.0000e+000 2.0000e+000 5.0213e+002 1.6000e+001 -1.9000e+000 1.0000e+000 2.0000e+000 5.0213e+002 1.6000e+001 改变初值x0后程序见Untitled1.m，输出结果：

-1.9000e+000 2.0000e+000 2.0000e+000 2.9923e+002 1.7000e+001 -1.9000e+000 2.0000e+000 2.0000e+000 2.9923e+002 1.7000e+001 -1.9000e+000 2.0000e+000 2.0000e+000 2.9923e+002 1.7000e+001 -1.9000e+000 2.0000e+000 2.0000e+000 2.9923e+002 1.7000e+001

-1.9000e+000 2.0000e+000 2.0000e+000 2.9923e+002 1.7000e+001 -1.9000e+000 2.0000e+000 2.0000e+000 2.9923e+002 1.7000e+001 -1.9000e+000 2.0000e+000 2.0000e+000 2.9923e+002 1.7000e+001 -1.9000e+000 2.0000e+000 2.0000e+000 2.9923e+002 1.7000e+001 结论：从2）4）输出观察和比较数据可以看出改变x0,最优解随之改变，而bfgs公式，混合2,3次插值，求得的局部最优解最接近全局最优解。

实验七T8：

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

3)minz(x12)2(x21)20.0425(x2x1) 12220.25x1x21解：3）在matlab中输入程序（exam07051grad.m，exam07051grad_run.m）：

function [f,g]=exam07051grad(x)

f=(x(1)-2)^2+(x(2)-1)^2+0.04/(-0.25*x(1)^2-x(2)^2+1)+5*(x(1)-2*x(2)+1)^2;% compute the function value at x

if nargout > 1 % fun called with 2 output arguments

g(1)=12*x(1) - 20*x(2) + x(1)/(50*(x(1)^2/4 + x(2)^2 - 1)^2) + 6 ; % compute the gradient evaluated at x

g(2)=42*x(2) - 20*x(1) + (2*x(2))./(25*(x(1)^2/4 + x(2)^2 - 1)^2) - 22; end

clear all; clc;

% comparing different algorithms: using gradient vector format short e x0=[-1.9,2];

'---case1: bfgs, hybrid 2,3 poly-------'

opt1=optimset('LargeScale','off', 'MaxFunEvals',1000,'GradObj','on'); [x1,v1,exit1,out1]=fminunc(@exam07051grad,x0,opt1) pause

'---case2: dfp, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','dfp');

[x2,v2,exit2,out2]=fminunc(@exam07051grad,x0,fopt) pause

'---case3: bfgs, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','bfgs');

[x3,v3,exit3,out3]=fminunc(@exam07051grad,x0,fopt) pause

'---case4: steep, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','steepdesc'); [x4,v4,exit4,out4]=fminunc(@exam07051grad,x0,fopt) pause

'---case5: bfgs, 3rd poly-------'

opt2=optimset(opt1,'LineSearchType','cubicpoly'); [x5,v5,exit5,out5]=fminunc(@exam07051grad,x0,opt2) pause

'---case6: dfp, 3rd poly-------' fopt=optimset(opt2,'HessUpdate','dfp');

[x6,v6,exit6,out6]=fminunc(@exam07051grad,x0,fopt) pause

'---case7:bfgs , 3rd poly-------' fopt=optimset(opt2,'HessUpdate','bfgs');

[x7,v7,exit7,out7]=fminunc(@exam07051grad,x0,fopt)

pause

'---case8: steep, 3rd poly-------'

fopt=optimset(opt2,'HessUpdate','steepdesc'); [x8,v8,exit8,out8]=fminunc(@exam07051grad,x0,fopt) pause

'++++ results of solutions ++++++' solutions=[x1;x2;x3;x4;x5;x6;x7;x8]; funvalues=[v1;v2;v3;v4;v5;v6;v7;v8];

iterations=[out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount;out7.funcCount;out8.funcCount]; [solutions,funvalues,iterations] 输出（太多，省略一部分）：

x8 =

1.7954e+000 1.3778e+000 v8 =

1.6904e-001

exit8 = 1

out8 =

iterations: 42

funcCount: 105

stepsize: 2.0518e-002 firstorderopt: 9.4180e-005

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x438 char] ans =

++++ results of solutions ++++++ ans =

1.7954e+000 1.3779e+000 1.6904e-001 9.0000e+000 1.7954e+000 1.3779e+000 1.6904e-001 1.1000e+001 1.7954e+000 1.3779e+000 1.6904e-001 9.0000e+000 1.7954e+000 1.3778e+000 1.6904e-001 1.0500e+002 1.7954e+000 1.3779e+000 1.6904e-001 9.0000e+000 1.7954e+000 1.3779e+000 1.6904e-001 1.1000e+001 1.7954e+000 1.3779e+000 1.6904e-001 9.0000e+000 1.7954e+000 1.3778e+000 1.6904e-001 1.0500e+002

clear all; clc;

% comparing different algorithms: using gradient vector format short e x0=[-1,2];

'---case1: bfgs, hybrid 2,3 poly-------'

opt1=optimset('LargeScale','off', 'MaxFunEvals',1000,'GradObj','on');

[x1,v1,exit1,out1]=fminunc(@exam07051grad,x0,opt1) pause

'---case2: dfp, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','dfp');

[x2,v2,exit2,out2]=fminunc(@exam07051grad,x0,fopt) pause

'---case3: bfgs, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','bfgs');

[x3,v3,exit3,out3]=fminunc(@exam07051grad,x0,fopt) pause

'---case4: steep, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','steepdesc'); [x4,v4,exit4,out4]=fminunc(@exam07051grad,x0,fopt) pause

'---case5: bfgs, 3rd poly-------'

opt2=optimset(opt1,'LineSearchType','cubicpoly'); [x5,v5,exit5,out5]=fminunc(@exam07051grad,x0,opt2) pause

'---case6: dfp, 3rd poly-------' fopt=optimset(opt2,'HessUpdate','dfp');

[x6,v6,exit6,out6]=fminunc(@exam07051grad,x0,fopt) pause

'---case7:bfgs , 3rd poly-------' fopt=optimset(opt2,'HessUpdate','bfgs');

[x7,v7,exit7,out7]=fminunc(@exam07051grad,x0,fopt) pause

'---case8: steep, 3rd poly-------'

fopt=optimset(opt2,'HessUpdate','steepdesc'); [x8,v8,exit8,out8]=fminunc(@exam07051grad,x0,fopt) pause

'++++ results of solutions ++++++' solutions=[x1;x2;x3;x4;x5;x6;x7;x8]; funvalues=[v1;v2;v3;v4;v5;v6;v7;v8];

iterations=[out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount;out7.funcCount;out8.funcCount]; [solutions,funvalues,iterations]

改变x0输入程序如下（ljxx0.m）：

clear all; clc;

% comparing different algorithms: using gradient vector format short e x0=[-1,2];

'---case1: bfgs, hybrid 2,3 poly-------'

opt1=optimset('LargeScale','off', 'MaxFunEvals',1000,'GradObj','on'); [x1,v1,exit1,out1]=fminunc(@exam07051grad,x0,opt1) pause

'---case2: dfp, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','dfp');

[x2,v2,exit2,out2]=fminunc(@exam07051grad,x0,fopt) pause

'---case3: bfgs, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','bfgs');

[x3,v3,exit3,out3]=fminunc(@exam07051grad,x0,fopt) pause

'---case4: steep, hybrid 2,3 poly-------' fopt=optimset(opt1,'HessUpdate','steepdesc');

[x4,v4,exit4,out4]=fminunc(@exam07051grad,x0,fopt) pause

'---case5: bfgs, 3rd poly-------'

opt2=optimset(opt1,'LineSearchType','cubicpoly'); [x5,v5,exit5,out5]=fminunc(@exam07051grad,x0,opt2) pause

'---case6: dfp, 3rd poly-------' fopt=optimset(opt2,'HessUpdate','dfp');

[x6,v6,exit6,out6]=fminunc(@exam07051grad,x0,fopt) pause

'---case7:bfgs , 3rd poly-------' fopt=optimset(opt2,'HessUpdate','bfgs');

[x7,v7,exit7,out7]=fminunc(@exam07051grad,x0,fopt) pause

'---case8: steep, 3rd poly-------'

fopt=optimset(opt2,'HessUpdate','steepdesc'); [x8,v8,exit8,out8]=fminunc(@exam07051grad,x0,fopt) pause

'++++ results of solutions ++++++' solutions=[x1;x2;x3;x4;x5;x6;x7;x8]; funvalues=[v1;v2;v3;v4;v5;v6;v7;v8];

iterations=[out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount;out7.funcCount;out8.funcCount]; [solutions,funvalues,iterations]

因太过繁琐，只输出最后x的值，输出： x8 =

7.8606e-001 8.2173e-001 v8 =

1.8420e+000

exit8 = 5

out8 =

iterations: 8 funcCount: 111

stepsize: 1.2680e-002 firstorderopt: 9.4132e-001

algorithm: 'medium-scale: Quasi-Newton line search' message: [1x362 char] ans =

++++ results of solutions ++++++ ans =

6.3105e-001 9.4892e-001 -1.8014e+014 1.5400e+002 1.7954e+000 1.3779e+000 1.6904e-001 2.5000e+001 6.3105e-001 9.4892e-001 -1.8014e+014 1.5400e+002 7.8606e-001 8.2173e-001 1.8420e+000 1.1100e+002 6.3105e-001 9.4892e-001 -1.8014e+014 1.5400e+002 1.7954e+000 1.3779e+000 1.6904e-001 2.5000e+001 6.3105e-001 9.4892e-001 -1.8014e+014 1.5400e+002 7.8606e-001 8.2173e-001 1.8420e+000 1.1100e+002

结论：从输出观察和比较数据可以看出改变x0,最优解随之改变，而bfgs公式，混合2,3次插值，求得的局部最优解最接近全局最优解。

实验七T5:

某分子由25 个原子组成，并且已经通过实验测量得到了其中某些原子对之间的距离（假设在平面结构上讨论），如下表所示。请你确定每个原子的位置关系。

原子对 距离 原子对 距离 原子对 距离 原子对 距离 (4,1) 0.9607 (5,4) 0.4758 (18,8) 0.8363 (15,13) 0.5725 (12,1) 0.4399 (12,4) 1.3402 (13,9) 0.3208 (19,13) 0.7660 (13,1) 0.8143 (24,4) 0.7006 (15,9) 0.1574 (15,14) 0.4394 (17,1) 1.3765 (8,6) 0.4945 (22,9) 1.2736 (16,14) 1.0952 (21,1) 1.2722 (13,6) 1.0559 (11,10) 0.5781 (20,16) 1.0422 (5,2) 0.5294 (19,6) 0.6810 (13,10) 0.9254 (23,16) 1.8255 (16,2) 0.6144 (25,6) 0.3587 (19,10) 0.6401 (18,17) 1.4325 (17,2) 0.3766 (8,7) 0.3351 (20,10) 0.2467 (19,17) 1.0851 (25,2) 0.6893 (14,7) 0.2878 (22,10) 0.4727 (20,19) 0.4995 (5,3) 0.9488 (16,7) 1.1346 (18,11) 1.3840 (23,19) 1.2277 (20,3) 0.8000 (20,7) 0.3870 (25,11) 0.4366 (24,19) 1.1271 (21,3) 1.1090 (21,7) 0.7511 (15,12) 1.0307 (23,21) 0.7060 (24,3) 1.1432 (14,8) 0.4439 (17,12) 1.3904 (23,22) 0.8052

解：

（1）建立模型

因为问题是要求出各原子间的距离关系，因此可以利用平面坐标轴来表示各个原

,则选择第1个点为参考点0,0󰀀子的相对位置，设第i个点的坐标为xi,yi 󰀀，所

以这些点的距离应该满足

Zmin[(xixj)2(yiyj)2dij2]2

i,j其中,dij是i原子和j原子之间的距离，此时，问题转化为无约束优化：

Zmin[(xixj)2(yiyj)2dij2]2

i,j（2）求解思路：

原题目中的问题转化为求无约束优化问题，可以借助matlab软件中的fminunc命令对Zmin[(xixj)2(yiyj)2dij2]2求出最优解。

i,j编程如下：

function f=juli1fun(x,d)

f(1)=(x(1,3))^2+(x(2,3))^2-d(1)^2; 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,1)-x(1,4))^2+(x(2,1)-x(2,4))^2-d(6)^2; f(7)=(x(1,1)-x(1,15))^2+(x(2,1)-x(2,15))^2-d(7)^2; f(8)=(x(1,1)-x(1,16))^2+(x(2,1)-x(2,16))^2-d(8)^2; f(9)=(x(1,1)-x(1,24))^2+(x(2,1)-x(2,24))^2-d(9)^2; f(10)=(x(1,2)-x(1,4))^2+(x(2,2)-x(2,4))^2-d(10)^2; f(11)=(x(1,2)-x(1,19))^2+(x(2,2)-x(2,19))^2-d(11)^2; f(12)=(x(1,2)-x(1,20))^2+(x(2,2)-x(2,20))^2-d(12)^2; f(13)=(x(1,2)-x(1,23))^2+(x(2,2)-x(2,23))^2-d(13)^2; f(14)=(x(1,3)-x(1,4))^2+(x(2,3)-x(2,4))^2-d(14)^2; f(15)=(x(1,3)-x(1,11))^2+(x(2,3)-x(2,11))^2-d(15)^2;

f(16)=(x(1,3)-x(1,23))^2+(x(2,3)-x(2,23))^2-d(16)^2; f(17)=(x(1,5)-x(1,7))^2+(x(2,5)-x(2,7))^2-d(17)^2; f(18)=(x(1,5)-x(1,12))^2+(x(2,5)-x(2,12))^2-d(18)^2; f(19)=(x(1,5)-x(1,18))^2+(x(2,5)-x(2,18))^2-d(19)^2; f(20)=(x(1,5)-x(1,24))^2+(x(2,5)-x(2,24))^2-d(20)^2; f(21)=(x(1,6)-x(1,7))^2+(x(2,6)-x(2,7))^2-d(21)^2; f(22)=(x(1,6)-x(1,13))^2+(x(2,6)-x(2,13))^2-d(22)^2; f(23)=(x(1,6)-x(1,15))^2+(x(2,6)-x(2,15))^2-d(23)^2; f(24)=(x(1,6)-x(1,19))^2+(x(2,6)-x(2,19))^2-d(24)^2; f(25)=(x(1,6)-x(1,20))^2+(x(2,6)-x(2,20))^2-d(25)^2; f(26)=(x(1,7)-x(1,13))^2+(x(2,7)-x(2,13))^2-d(26)^2; f(27)=(x(1,7)-x(1,17))^2+(x(2,7)-x(2,17))^2-d(27)^2; f(28)=(x(1,8)-x(1,12))^2+(x(2,8)-x(2,12))^2-d(28)^2; f(29)=(x(1,8)-x(1,14))^2+(x(2,8)-x(2,14))^2-d(29)^2; f(30)=(x(1,8)-x(1,21))^2+(x(2,8)-x(2,21))^2-d(30)^2; f(31)=(x(1,9)-x(1,10))^2+(x(2,9)-x(2,10))^2-d(31)^2; f(32)=(x(1,9)-x(1,12))^2+(x(2,9)-x(2,12))^2-d(32)^2; f(33)=(x(1,9)-x(1,18))^2+(x(2,9)-x(2,18))^2-d(33)^2; f(34)=(x(1,9)-x(1,19))^2+(x(2,9)-x(2,19))^2-d(34)^2; f(35)=(x(1,9)-x(1,21))^2+(x(2,9)-x(2,21))^2-d(35)^2; f(36)=(x(1,10)-x(1,17))^2+(x(2,10)-x(2,17))^2-d(36)^2; f(37)=(x(1,10)-x(1,24))^2+(x(2,10)-x(2,24))^2-d(37)^2; f(38)=(x(1,11)-x(1,14))^2+(x(2,11)-x(2,14))^2-d(38)^2; f(39)=(x(1,11)-x(1,16))^2+(x(2,11)-x(2,16))^2-d(39)^2; f(40)=(x(1,12)-x(1,14))^2+(x(2,12)-x(2,14))^2-d(40)^2; f(41)=(x(1,12)-x(1,18))^2+(x(2,12)-x(2,18))^2-d(41)^2; f(42)=(x(1,13)-x(1,14))^2+(x(2,13)-x(2,14))^2-d(42)^2; f(43)=(x(1,13)-x(1,15))^2+(x(2,13)-x(2,15))^2-d(43)^2; f(44)=(x(1,15)-x(1,19))^2+(x(2,15)-x(2,19))^2-d(44)^2; f(45)=(x(1,15)-x(1,22))^2+(x(2,15)-x(2,22))^2-d(45)^2;

f(46)=(x(1,16)-x(1,17))^2+(x(2,16)-x(2,17))^2-d(46)^2; f(47)=(x(1,16)-x(1,18))^2+(x(2,16)-x(2,18))^2-d(47)^2; f(48)=(x(1,18)-x(1,19))^2+(x(2,18)-x(2,19))^2-d(48)^2; f(49)=(x(1,18)-x(1,22))^2+(x(2,18)-x(2,22))^2-d(49)^2; f(50)=(x(1,18)-x(1,23))^2+(x(2,18)-x(2,23))^2-d(50)^2; f(51)=(x(1,20)-x(1,22))^2+(x(2,20)-x(2,22))^2-d(51)^2; f(52)=(x(1,21)-x(1,22))^2+(x(2,21)-x(2,22))^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.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]'; x0=[zeros(1,24);ones(1,24)];

[x,norms]=lsqnonlin(@juli1fun,x0,[],[],[],d) p=x;

a=[0,x(1,:)]; b=[0,x(2,:)]; plot(a,b,'*')

则输出的结果是 x =

Columns 1 through 6

0.7497 1.8726 0.9484 0.9462 1.4904 0.7934 -1.0399 -0.7569 -0.0931 -0.5578 -0.1510 -0.3849

Columns 7 through 12

1.1476 0.4874 1.3734 -0.5282 -0.2585 -0.1734

Columns 13 through 18

0.6863 0.6291 1.2468 -0.4824 -0.0522 -1.4191

Columns 19 through 24

1.1737 1.2655 1.7338 -0.3692 0.1777 -0.5092 norms =

0.0307 图像：

0.7966 -0.1222 0.3911 -1.3192 1.9656 0.2582 -0.3925 0.5555 -0.1711 -0.6129 1.8013 1.3289 -1.0690 -0.7964 1.5060 1.1416 0.3220 -0.4337

（3）数据分析：

由上述数据可分别得出原子的坐标位置，和原子平面分布的状况。根据拟合的范数norms>0,说明拟合的结果可靠。 （4）引出其他问题：

上述程序中，初值是随意给定的，但改变初值时，结果会发生改变。

实验七T7：

经济学中著名的柯布-道格拉斯（Cobb-Douglas）生产函数的一般形式为

Q(K,L)aKL,0,1

其中Q, K, L分别表示产值、资金、劳动力，式中,,a要由经济统计数据确定.现有《中国统计年鉴（2003）》给出的统计数据如表所示，请分别用线性和非线性最小二乘法拟合求出式中,,a的并解释,的含义。 年份 1984 1985 总产值/万亿元 0.7171 0.8964 资金/万亿元 0.0910 0.2543 劳动力/亿人 4.8179 4.9873 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 1.0202 1.1962 1.4928 1.6909 1.8548 2.1618 2.6638 3.4634 4.6759 5.8478 6.7885 7.4463 7.8345 8.2068 9.9468 9.7315 10.4791 0.3121 0.3792 0.4754 0.4410 0.4517 0.5595 0.8080 1.3072 1.7042 2.0019 2.2914 2.4941 2.8406 2.9854 3.2918 3.7314 4.3500 5.1282 5.2783 5.4334 5.5329 6.4749 6.5491 6.6152 6.6808 6.7455 6.8065 6.8950 6.9820 7.0637 7.1394 7.2085 7.3025 7.3740 其中总产值取自“国内生产总值”，资金取自“固定资产投资”，劳动力取自“就业人员” 解：

先进行模型分析： 题目已经确定了此数据的拟合形式为Cobb-Douglas生产函数方程，再以这些数据分别给出线性与非线性的最小二乘拟合，则可以确定式中的α,β,a。 在构造函数时，以x向量的分量分别表示a,α,β,以t矩阵的第一列表示资金K，第二列表示劳动力L。

利用Matlab编写非线性拟合程序代码，具体程序如下： function f=cobb(x,t)

f=x(1).*(t(1,:).^x(2)).*(t(2,:).^x(3)); end

x0=[1 1 1];

K=[0.0910 0.2543 0.3121 0.3792 0.4754 0.4410 0.4517 0.5595 0.8080 1.3072 1.7042 2.0019

2.2914 2.4941 2.8406 2.9854 3.2918 3.7314 4.3500];

L=[4.8179 4.9873 5.1282 5.2783 5.4334 5.5329 6.4749 6.5491 6.6152 6.6808 6.7455 6.8065

6.8950 6.9820 7.0637 7.1394 7.2085 7.3025 7.3740];

Q=[0.7171 0.8964 1.0202 1.1962 1.4928 1.6909 1.8548 2.1618 2.6638 3.4634 4.6759 5.8478

6.7885 7.4463 7.8345 8.2068 9.9468 9.7315 10.4791];

[x,norm,res,ef,out]=lsqcurvefit(@cobb,x0,[K;L],Q,[0,0,0],[1,1,1])

输出结果如下： x =

0.83444 0.7736 0.73117 norm = 2.7489 res =

Columns 1 through 4

-0.30483 0.040363 0.099976 0.13387 Columns 5 through 8

0.12539 -0.14366 -0.086582 -0.057732 Columns 9 through 12

0.15279 0.65269 0.41327 -0.045543 Columns 13 through 16

-0.28605 -0.43922 -0.019363 -0.021674 Columns 17 through 19

-1.0567 0.15712 0.73501

ef = 3 out =

firstorderopt: 1.7046e-005 iterations: 5 funcCount: 24 cgiterations: 0

algorithm: 'trust-region-reflective' message: [1x460 char]

非线性拟合给出的结果，a=0.83444,α=0.7736,β=0.73117，其误差平方和norm=2.7489，稍大。在接下来的分析中将看到它与原数据的拟合程度。 再编写相关程序作图观察QKL的趋势： t=1984:2002; plot(t,Q,'r'); gtext('Q');

hold on,plot(t,K,'k'),hold off gtext('K');

hold on,plot(t,L,'b'),hold off gtext('L')

从形式上观察，我们不妨将生产函数适当对数处理，将所给的Cobb-Douglas方程两边取对数，则可将非线性方程化为线性形式：

lnQ(K,L)lnalnKlnL再令ylnQ,klnK,llnL,mlna，则方

程化为线性形式：

ymkl

令

则问题转化为求超定方程：xy根据最小二乘准则计算可以得到矩阵中的待定系数。由于要求线性拟合，故采用lsqlin函数编写如下线性拟合代码： K=[0.0910 0.2543 0.3121 0.3792 0.4754 0.4410 0.4517 0.5595 0.8080 1.3072 1.7042 2.0019 2.2914 2.4941 2.8406 2.9854 3.2918 3.7314 4.3500]';

L=[4.8147 4.9873 5.1282 5.2783 5.4334 5.5329 6.4749 6.5491 6.6152 6.6808 6.7455 6.8065 6.8950 6.9820 7.0637 7.1394 7.2085 7.3025 7.3740]';

Q=[0.7171 0.8964 1.0202 1.1962 1.4928 1.6909 1.8548 2.1618 2.6638 3.4634 4.6759 5.8478 6.7885 7.4463 7.8345 8.2068 9.9468 9.7315 10.4791]';

k=log(K); l=log(L); f=log(Q);

M=[ ones(19,1),k,l];

[x,resnorm,residual,exit,out]=lsqlin(M,f,[],[],[],[],[-5,0,0]',[0,1,1]') a=exp(x(1)) alpha=x(2) beta=x(3)

输出结果如下： x =

-0.0066453 0.73489 0.64329 resnorm = 0.3578 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

直接左除 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)

输出结果如下： a =

0.31777 alpha = 0.66054 beta =

1.2627

不符合实际情况，应是左除无法加入参数条件限制导致的。下对此情况不予讨论。 作图对比 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('原始数据与拟合函数对比')

黑线为原始数据的折线，红色为非线性拟合，绿色为取对数后的线性拟合。由于原始数据只有19个，量比较小，还不足以看出拟合的效果究竟如何。从图上看，

它们都基本描述出了Q的走势，理论上讲，线性拟合的效果更好一些。 最后结果分析： 方法 非线性拟合 对数线性拟合 a 0.83444 0.9938 α 0.7736 0.73489 β norm iter 5 6 0.73117 2.4789 0.64329 0.35783 从数值上看，二者计算出的结果相差比较大。但是从图上看没有太大问题。数据分析，对数线性拟合的norm=0.3579，比起用非线性拟合的数据norm=2.4789，此norm值要小得多，这说明用取对数后用线性拟合效果很好。而且在试探的过程中，发现结果对给定参数的范围的变化较敏感。但是对初值的变化很不敏感。

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

数学实验作业第五组 侯永利（121170011） 吕菊香（121170012） 陈智杰（121170018） 2014年8月10日