數(shù)學(xué)實(shí)驗(yàn)8月13日作業(yè)

1、.1.取不同的初值計算下列平方和形式的非線性規(guī)劃，盡可能求出所有局部極小點(diǎn)，進(jìn)而找出全局極小點(diǎn)，并對不同算法（搜索方向、搜索步長、數(shù)值梯度與分析梯度等）的結(jié)果進(jìn)行分析、比較。(2). ,(4).，其中.解：(2)由題意設(shè)先求（文件名為a.m），再求（文件名為b.m）,解得近似解為或或或或或或或或,取初值（exam0705grad_run.m）利用matlab計算可得以下表中數(shù)據(jù)：搜索方向步長搜索最優(yōu)解最優(yōu)解最優(yōu)值目標(biāo)函數(shù)調(diào)用次數(shù)Case 1bfgs混合二、三次插值-2.8996e+0002.1584e-0014.2939e-0065.0000e+000Case 2dfp-2.8996e+000

2、2.1584e-0014.2939e-0065.0000e+000Case 3bfgs-2.8996e+0002.1584e-0014.2939e-0065.0000e+000Case 4steepdesc-2.8996e+0002.1584e-0014.2939e-0065.0000e+000Case 5bfgs三次插值-2.8996e+0002.1584e-0014.2939e-0065.0000e+000Case 6dfp-2.8996e+0002.1584e-0014.2939e-0065.0000e+000Case 7bfgs-2.8996e+0002.1584e-0014.2939

3、e-0065.0000e+000Case 8steepdesc-2.8996e+0002.1584e-0014.2939e-0065.0000e+000當(dāng)取初值（exam07051grad_run.m）利用matlab計算可得以下表中數(shù)據(jù)：搜索方向步長搜索最優(yōu)解最優(yōu)解最優(yōu)值目標(biāo)函數(shù)調(diào)用次數(shù)Case 1bfgs混合二、三次插值3.4578e+000-7.9722e-0028.3883e-0121.4000e+001Case 2dfp3.4578e+000-7.9722e-0027.2387e-0101.6000e+001Case 3bfgs3.4578e+000-7.9722e-0028.388

4、3e-0121.4000e+001Case 4steepdesc3.4578e+000-7.9722e-0025.0577e-009 1.7000e+001Case 5bfgs三次插值3.4578e+000-7.9722e-0028.3883e-0121.4000e+001Case 6dfp3.4578e+000-7.9722e-0027.2387e-0101.6000e+001Case 7bfgs3.4578e+000-7.9722e-0028.3883e-0121.4000e+001Case 8steepdesc3.4578e+000-7.9722e-0025.0577e-009 1.70

5、00e+001因?yàn)椴煌介L搜索到相同的最優(yōu)解，但最優(yōu)值不同，所以此時初值離最優(yōu)解很遠(yuǎn)。繼續(xù)取不同初值求解（只需在程序中不斷改變初值即可）。當(dāng)初值取時可得到以下數(shù)據(jù)搜索方向步長搜索最優(yōu)解最優(yōu)解最優(yōu)值目標(biāo)函數(shù)調(diào)用次數(shù)Case 1bfgs混合二、三次插值-2.1181e+001-3.6468e+0011.3691e-0061.7000e+001Case 2dfp-2.1158e+001-3.6504e+0011.9079e+0001.9500e+002Case 3bfgs-2.1181e+001-3.6468e+0011.3691e-0061.7000e+001Case 4steepdesc-2.1

6、174e+001-3.6478e+0011.5885e-001 7.7100e+002Case 5bfgs三次插值-2.1181e+001-3.6468e+0011.3691e-0061.7000e+001Case 6dfp-2.1158e+001-3.6504e+0011.9079e+0001.9500e+002Case 7bfgs-2.1181e+001-3.6468e+0011.3691e-0061.7000e+001Case 8steepdesc-2.1174e+001-3.6478e+0011.5885e-001 7.7100e+002綜上全局極小點(diǎn)為，局部極小點(diǎn)，。（2）8 取不同

7、的初值計算下列非線性規(guī)劃，盡可能求出所有局部極小點(diǎn)，進(jìn)而找出全局極小點(diǎn)，并對不同算法（搜索方向、搜索步長、數(shù)值梯度與分析梯度等）的結(jié)果進(jìn)行分析、比較。 (2),其中，,。解：(2)1. 設(shè)計程序如exam070702fun.m.建立功能函數(shù)，便于調(diào)用。functionf,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 nargout1 g(1)=2*(x(1)-a1(1)/(x(1)-a1(1)2+(x(2)-a1(2)2+c(1)2+2*(x(1)-a2

8、(1)/(x(1)-a2(1)2+(x(2)-a2(2)2+c(2)2 %該為梯度函數(shù)的兩個方向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)2end設(shè)計腳本程序如exam070702.m，opt1=optimset(LargeScale,off,MaxfunEvals,1000,Tolfun,1e-8,TolX,1e-8,Gradobj,on); opt2=optimset(opt1,HessUpdate,dfp); opt3=optimset(op

9、t1,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

10、,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)運(yùn)行結(jié)果及分析。得出答案；制成表格初值設(shè)定為（20，43）搜索方向 步長搜索最優(yōu)解最優(yōu)值 目標(biāo)函數(shù)調(diào)用次數(shù)BFGS混合二、三次插值0.2024e-004

11、 0.4775e-0042.0423e-00512DFP混合二、三次插值0.2024e-004 0.4775e-0042.0423e-00512最速下降混合二、三次插值0.2024e-004 0.4775e-0042.0423e-00512BFGS三次插值0.2024e-004 0.4775e-0042.0423e-00512DFP三次插值0.2024e-004 0.4775e-0042.0423e-00512最速下降三次插值0.2024e-004 0.4775e-0042.0423e-00512 可以看到采用分析方法(給定梯度函數(shù)) 比起默認(rèn)的數(shù)值方法得到結(jié)果所需要的調(diào)用次數(shù)明顯減少。DFP

12、算法的函數(shù)調(diào)用次數(shù)較少。 改變初值，采用默認(rèn)搜索方向BFGS，默認(rèn)搜索步長，給定梯度函數(shù)，得到的結(jié)果如下:X(1)X(2)最優(yōu)值迭代次數(shù)(3,3)2.49333.7345-1.03292(4,4)3.67113.95610.76401(5,5)5 50.25161(2,7)2.98994.1982-0.31781可以看到，改變初值對結(jié)果的影響很大.結(jié)論： 由上面計算結(jié)果得出，無論采用那種計算方法最終都可以得到函數(shù)的極小值的，而采用DFP算法函數(shù)的調(diào)用次數(shù)較少；分析法計算梯度函數(shù)效果較好；最終結(jié)果取哪一個局部極值同初值密切相關(guān)。 5. 某分子由25 個原子組成，并且已經(jīng)通過實(shí)驗(yàn)測量得到了其中某些

13、原子對之間的距離（假設(shè)在平面結(jié)構(gòu)上討論），如下表所示。請你確定每個原子的位置關(guān)系。原子對距離原子對距離原子對距離原子對距離（4,1）0.9607(5,4)0.4758(18,8)0.8363(15,13)0.5725(12,1)0.4399(12,41.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.0

14、422(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

15、(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解:問題分析每個原子的位置都是未知的，在坐標(biāo)系中只有相對的位置參數(shù)，不妨固定原子1的坐標(biāo)為。并且分子可以在平面內(nèi)任意旋轉(zhuǎn)，很難確定每個原子的絕對位置，表中給出了52組數(shù)據(jù)，可以得到52個方程，而未知數(shù)個數(shù)為48個（每個原子的，坐標(biāo)）。屬于超定方程組，沒有確定的解，只能求得最優(yōu)解。故采用最小方差標(biāo)準(zhǔn)來求最小二乘解得到最優(yōu)解。模型建立設(shè)第個點(diǎn)所在的位置為，因?yàn)樗蟮氖歉髟娱g的位置關(guān)系，可以設(shè)定第一個

16、點(diǎn)坐標(biāo)為，然后再計算其他原子的位置，使它在最大程度上滿足上表中提供的數(shù)據(jù)，即讓達(dá)到最小，求得的坐標(biāo)值就是最優(yōu)的解。其中表示第個原子和第個原子之間的距離，數(shù)據(jù)如上表所示，問題轉(zhuǎn)化為無約束優(yōu)化；算法設(shè)計根據(jù)題目中所給的數(shù)據(jù)，利用matlab編寫函數(shù)文件（文件名為diatance.m）,再主程序?yàn)閑xam5.m,輸出結(jié)果每個原子的位置如下（即第二個原子的位置為（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.6

17、158 -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為了形象直觀地表示出各點(diǎn)的位置，畫

18、出如下的散點(diǎn)圖7. 經(jīng)濟(jì)學(xué)中著名的柯布道格拉斯（CobbDouglas）生產(chǎn)函數(shù)的一般形式為，其中,分別表示產(chǎn)值、資金、勞動力，式中,要由經(jīng)濟(jì)統(tǒng)計數(shù)據(jù)確定?，F(xiàn)有中國統(tǒng)計年鑒（2003）給出的統(tǒng)計數(shù)據(jù)如下表，請用非線性最小二乘擬合求出式中的，并解釋，的含義。年份總產(chǎn)值（萬億元）資金（萬億元）勞動力（億人）19840.71710.09104.817919850.89640.25434.987319861.02020.31215.128219871.19620.37925.278319881.49280.47545.433419891.69090.44105.532919901.85480.4517

19、6.474919912.16180.55956.549119922.66380.80806.615219933.46341.30726.680819944.67591.70426.745519955.84782.00196.806519966.78852.29146.895019977.44632.49416.982019987.83452.84067.063719998.20682.98547.139420009.94683.29187.208520019.73153.73147.3025200210.47914.35007.3740其中總產(chǎn)值取自“國內(nèi)生產(chǎn)總值”，資金取自“固定資產(chǎn)投資”，

20、勞動力取自“就業(yè)人員”。解：模型分析題目已經(jīng)確定了此數(shù)據(jù)的擬合形式為Cobb-Douglas生產(chǎn)函數(shù)方程，再以這些數(shù)據(jù)分別給出線性與非線性的最小二乘擬合，則確定式中的,。 在構(gòu)造函數(shù)時，以向量的分量分別表示,以矩陣的第一列表示資金，第二列表示勞動力。非線性擬合程序代碼先用mutlab編寫函數(shù)（文件名為examfun7.m）,主程序?yàn)閑xam7.m,輸出結(jié)果如下x = 0.8345 0.7736 0.7312norm = 2.7489res = Columns 1 through 8 -0.3048 0.0404 0.1000 0.1339 0.1254 -0.1437 -0.0866 -0.0

21、577 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.7350ef = 3out = firstorderopt: 4.7566e-004 iterations: 14 funcCount: 60 cgiterations: 0 algorithm: large-scale: trust-region reflective Newton message: 1x460 char非線性擬合給出的結(jié)果，,

22、,，其誤差平方和，稍大。在接下來的分析中將看到它與原數(shù)據(jù)的擬合程度。 線性擬合程序代碼1. 作圖觀察QKL的趨勢用mutlab編寫程序（文件名為exam71.m）,運(yùn)行，得到如下圖像2. 對數(shù)處理從形式上觀察，我們不妨將生產(chǎn)函數(shù)適當(dāng)處理，將所給的Cobb-Douglas方程兩邊取對數(shù)，則可將非線性方程化為線性形式：再令,則方程可以轉(zhuǎn)化為。令，則問題轉(zhuǎn)化為超定方程根據(jù)最小二乘準(zhǔn)則計算可以得到矩陣中的待定系數(shù)。由于要求線性擬合，故采用lsqlin 函數(shù)。 3. 線性擬合代碼用matlab編寫線性擬合代碼（文件名為exam72.m），輸出結(jié)果如下，x = -0.0066453 0.73489 0.6

23、4329 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.035

24、332 -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=Mf a=exp(s(1) alpha=s(2) beta=s(3)輸出結(jié)果為：s = -1.1464 0.6605 1.2627a = 0.3178alpha = 0.6605beta =1.2627不符合實(shí)際情況，應(yīng)是左除無法加入?yún)?shù)條件限制導(dǎo)致的。下對此情況不予討論4. 作圖對比程序?yàn)閑xam74.m,代碼如下：for t=1