数学试验8月13日作业

1、1.取不同的初值计算下列平方和形式的非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法（搜索方向、搜索步长、数值梯度与分析梯度等）的结果进行分析、比较。2(2).minx112x221149x1249x2281x122324x2681(4).min100x310X1,X22Xi2X2x32,其中1arctgXi,X2,1.-arctgX2.X1X2.X1解：（2）由题意设fx1,x2先求-,（文件名为）x1x2,X11一，Xi22Xi12X221149xi2-2c.49x281x12324x22681X1X23.46-或0.08X1X220.62或34.53X1X224.1

2、3或17.92X1X20.83或23.72X1X2XifX2X1X2（文件名为）0，解得近似解为X1X22.90_或0.2221.1836.470.17,取初值47.72X1X2X）0.280.282.90,0.22X1X222.37或17.49利用matlab计算可得以下表中数据:搜索方向步长搜索最优解X1最优解X2最优值目标函数调用次数Case1bfgs混合二、三次插值-21.181,-36+000.468+000Case2dfp+000+000Case3bfgs+000+000Case4steepde+000+000scCase5bfgs三次插值+000+000Case6dfp+000+

3、000Case7bfgs+000+000Case8steepdesc+000+000当取初值xo3.46,0.08()利用matlab计算可得以下表中数据:搜索方向步长搜索最优解X1最优解X2最优值目标函数调用次数Case1bfgs混合二、三次插值+000+001Case2dfp+000+001Case3bfgs+000+001Case4steepdesc+000+001Case5bfgs三次插值+000+001Case6dfp+000+001Case7bfgs+000+001Case8steepdesc+000+001因为不同步长搜索到相同的最优解，但最优值不同，所以此时初值离最优解很远。继

4、续取不同初值求解(只需在程序中不断改变初值即可)。当初值取x021.18,36.47时可得到以下数据搜索方步长搜最优解X1最优解X2最优值目标函数调向索用次数Case1bfgs混合二、三次插值+001+001+001Case2dfp+001+001+000+002Case3bfgs+001+001+001Case4steepdesc+001+001+002Case5bfgs三次插值+001+001+001Case6dfp+001+001+000+002Case7bfgs+001+001+001Case8steepdesc+001+001+002综上全局极小点为-21.181,-36.468，局

5、部极小点-2.8996,0.21584,-21.181,-36.468。8.取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法(搜索方向、搜索步长、数值梯度与分析梯度等)的结果进行分析、比较。/O112(2)minztT,x(xa1)(xa1)c1(xa2)(xa2)c2其中c=(0.7,0.73),a1=(4,4)T,a2=(2.5,3.8)T。解：1.设计程序如.建立功能函数，便于调用。functionf,g=exam070702fun(x)a1=4,4'a2=,'c=,;f=1/(x-a1)'*(x-a1)+c(1)-1/(

6、x-a2)'*(x-a2)+c(2);ifnargout>1该为梯度函数的两个方向g(1)=2*(x(1)-a1(1)/(x(1)-a1(1)A2+(x(2)-a1(2)A2+c(1)A2+2*(x(1)-a2)/(x(1)-a2(1)A2+(x(2)-a2(2)A2+c(2)A2%g(2)=2*(x(2)-a1(2)/(x(1)-a1(1)A2+(x(2)-a1(2)A2+c(1)A2+2*(x(2)-a2(2)/(x(1)-a2(1)A2+(x(2)-a2(2)A2+c(2)A2end设计脚本程序如,opt1=optimset('LargeScale',

7、9;off,'MaxfunEvals',1000,'Tolfun',1e-8,'ToIX',1e-8,'Gradobj','on');opt2=optimset(opt1,'HessUpdate','dfp');opt3=optimset(opt1,'HessUpdate','steepdesc');opt4=optimset(opt1,'lineSearchType','cubicpoly');opt5=optimse

8、t(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=fm

9、inunc(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)搜索方向步长搜索最优解最优值目标函数调用次数BFGS混合二、三次插值12DFP混合二、三次插值12最速卜降混合

10、二、三次插值12BFGS三次插值12DFP三次插值12最速卜降三次插值12可以看到采用分析方法（给定梯度函数）比起默认的数值方法得到结果所需要的调用次数明显减少。DFPJ法的函数调用次数较少。改变初值，采用默认搜索方向BFGS默认搜索步长，给定梯度函数，得到的结果如下：X(1)X(2)最优值迭代次数(3.3)2(4,4)1(5.5)551(2,7)1可以看到，改变初值对结果的影响很大.结论：由上面计算结果得出，无论采用那种计算方法最终都可以得到函数的极小值的，而采用DFP：法函数的调用次数较少；分析法计算梯度函数效果较好；最终结果取哪一个局部极值同初值密切相关。5.某分子由25个原子组成，并且

11、已经通过实验测量得到了其中某些原子对之间的距离（假设在平面结构上讨论），如下表所示。请你确定每个原子的位置关系原子对距离原距离原子距离原子距离子对对对(4,1)(5,4(18,8(15,1)3)(12,1)(12,(13,9(19,14)3)(13,1)(24,(15,9(15,14)4)(17,1)(8,6(22,9(16,1)4)(21,1)(13,(11,1(20,16)0)6)(5,2)(19,(13,1(23,16)0)6)(16,2)(25,(19,1(18,16)0)7)(17,2)(8,7(20,1(19,1)0)7)(25,2)(14,(22,1(20,17)0)9)(5,3

12、)(16,(18,1(23,17)1)9)(20,3)(20,(25,1(24,17)1)9)(21,3)(21,(15,1(23,27)2)1)(24,3)(14,8)(17,12)(23,22)解:问题分析每个原子的位置都是未知的，在坐标系中只有相对的位置参数，不妨固定原子1的坐标为(0,0)。并且分子可以在平面内任意旋转，很难确定每个原子的绝对位置，表中给出了52组数据，可以得到52个方程，而未知数个数为48个(每个原子的x,y坐标)。属于超定方程组，没有确定的解，只能求得最优解。故采用最小方差标准来求最小二乘解得到最优解。模型建立设第i个点所在的位置为(x，yi),因为所求的是各原子间

13、的位置关系，可以设定第一个点坐标为(Xi,yJ(0,0),然后再计算其他原子的位置(x,yj,使它在最大2程度上满足上表中提供的数据，即让z=(XiXi)2(yiy。)2dij2达到最小,i.j求得的坐标值就是最优的解。其中dj表示第i个原子和第j个原子之间的距离，数据如上表所示，问题转化为无约束优化；min=(xXi)2(yiy°)2dj2i.j算法设计根据题目中所给的数据，利用matlab编写函数文件(文件名为)，再主程序为，输出结果每个原子的位置如下(即第二个原子的位置为(，)，第三个原子为(1,6142,1,1211),)为了形象直观地表示出各点的位置，画出如下的散点图-0.

14、2-0.4-0.6-0.8-1-1.2-1.4-1.6-1.8-0.50.51.57.经济学中著名的柯布一道格拉斯（Cobb-Douglas）生产函数的一般形式为Q（K,L）aKL,0,1其中Q,K,L分别表示产值、资金、劳动力，式中，a要由经济统计数据确定。现有中国统计年鉴（2003）»给出的统计数据如下表，请用非线性最小二乘拟合求出式中的，a,并解释，的含义。年份总产值（万亿元）资金（万亿元）力动力（亿人）1984198519861987198819891990199119921993199419951996199719981999200020012002其中总产值取自“国内生产

15、总值”，资金取自“固定资产投资”，劳动力取自“就业人员”。解：模型分析题目已经确定了此数据的拟合形式为Cobb-Douglas生产函数方程，再以这些数据分别给出线性与非线性的最小二乘拟合，则确定式中的a,。在构造函数时，以x向量的分量分别表示a,，以t矩阵的第一列表示资金K,第二列表示劳动力L。非线性拟合程序代码先用mutlab编写函数（文件名为）,主程序为，输出结果如下norm=res=Columns1through8Columns9through16Columns17through19ef=3out=firstorderopt:iterations:14funcCount:60cgiter

16、ations:0algorithm:'large-scale:trust-regionreflectiveNewtonmessage:1x460char非线性拟合给出的结果，a=0.83444,=0.7736,=0.73117,其误差平方和norm=2.7489,稍大。在接下来的分析中将看到它与原数据的拟合程度。线性拟合程序代码1 .作图观察QKL勺趋势用mutlab编写程序(文件名为)，运行，得到如下图像2 .对数处理从形式上观察，我们不妨将生产函数适当处理，将所给的Cobb-Douglas方程两边取对数，则可将非线性方程化为线性形式：lnQ(K,L)lnalnKlnL再令ylnQ,

17、llnLmlna,则方程可以转化为ymkl01k1令1k(a,)T,v(vi,V2,V3|,Vi9)1k19l19则问题转化为超定方程根据最小二乘准则计算可以得到矩阵中的待定系数。由于要求线性拟合，故采用Isqlin函数。3 .线性拟合代码用matlab编写线性拟合代码（文件名为），输出结果如下,x=resnorm=residual=Columns1through4Columns5through8Columns9through12Columns13through16Columns17through19exit=3out=iterations:6algorithm:1x43charfirstor

18、deropt:cgiterations:5message:1x123chara=alpha=beta=4.直接左除程序名为，代码如下：k=log(K);l=log(L);f=log(Q);M=ones(19,1),k,l;s=Mfa=exp(s(1)alpha=s(2)beta=s(3)输出结果为：s=alpha=beta=卜对此情况不予讨论不符合实际情况，应是左除无法加入参数条件限制导致的4.作图对比程序为，代码如下：fort=1:19Q1(t)=*K(t)A*L(t)A;endfori=1:19Q2(i)=*K(i)A*L(i)A；endts=1984:2002;plot(ts,Q,'k'),holdon原始数据与拟合函数对比,)plot(ts,Q1,'r'),holdonplot