MATLAB的黃金分割法與拋物線插值法

...wd......wd......wd...MATLAB的黃金分割法與拋物線插值法摘要為了求解最優(yōu)化模型的最優(yōu)解，可使用基于MATLAB算法編程的黃金分割法與拋物線插值法，來實現求解的過程。黃金分割法是通過所選試點的函數值而逐步縮短單谷區(qū)間來搜索最優(yōu)點，利用迭代進而得出結論。拋物線插值法亦稱二次插值法，是一種多項式插值法，逐次以擬合的二次曲線的極小點，逼近原尋求函數極小點的一種方法。通過將MATLAB與最優(yōu)化問題相結合，不僅可以加深對黃金分割法、拋物線插值法的根本理解和算法框圖及其步驟的全面理解，也有利于幫助我們掌握MATLAB的使用方法。關鍵詞：MATLAB，黃金分割法，拋物線插值法，最優(yōu)解，迭代目錄1.黃金分割法?????????????????????????????????????????????????????????????????21.1算法原理??????????????????????????????????????????????????????21.2算法步驟??????????????????????????????????????????????????????21.3黃金分割法算法框圖???????????????????????????????????????32.拋物線插值法??????????????????????????????????????????????????????4 2.1算法原理??????????????????????????????????????????????????????4 2.2算法步驟??????????????????????????????????????????????????????4 2.3拋物線插值法算法框圖?????????????????????????????????????53.算法的MATLAB實現???????????????????????????????????????????????63.1黃金分割法程序代碼??????????????????????????????????????????????63.2實例驗證??????????????????????????????????????????????????????63.3誤差分析??????????????????????????????????????????????????????93.4拋物線插值法程序代碼?????????????????????????????????????103.5實例驗證??????????????????????????????????????????????????????113.6誤差分析??????????????????????????????????????????????????????123.7黃金分割法與拋物線插值法的比照?????????????????????????13引言為了對最優(yōu)化問題進展求解，為了更好的掌握MATLAB的運用，本文主要介紹一維最優(yōu)化問題的一維搜索法黃金分割法與拋物線插值法，利用MATLAB算法將其實現。黃金分割法也叫0.618法，它是一種基于區(qū)間收縮的極小點搜索算法，當用進退法確定搜索區(qū)間后，我們只知道極小點包含于搜索區(qū)間內，但是具體是哪個點，無法得知。黃金分割法的思想很直接，既然極小點包含于搜索區(qū)間內，那么可以不斷的縮小搜索區(qū)間，就可以使搜索區(qū)間的端點逼近到極小點。拋物線插值法(parabolicinterpolationmet-hod)亦稱二次插值法一種多項式插值法.逐次以擬合的二次曲線的極小點，逼近原尋求函數極小點的一種方法,能求得精度較高的解，但速度會比照慢。時至今日，經過MathWorks公司的不斷完善，MATLAB已經開展成為適合多學科、多種工作平臺的功能強勁的大型軟件。在國外，MATLAB已經經受了多年考驗。在歐美等高校，MATLAB已經成為線性代數、自動控制理論、數理統計、數字信號處理、時間序列分析、動態(tài)系統仿真等高級課程的根本教學工具；成為攻讀學位的大學生、碩士生、博士生必須掌握的根本技能。在設計研究單位和工業(yè)部門，MATLAB被廣泛用于科學研究和解決各種具體問題。1.黃金分割法1.1算法原理黃金分割法的根本原理：黃金分割法又稱0.618法，它是通過不斷縮短搜索區(qū)間的長度來尋求一維函數的極小點。這種方法的原理是：在搜索區(qū)間[a,b]內按如下規(guī)則對稱地取兩點：計算它們的函數值f1=f(a1),f2=f(a2),比照它們的大小，結果有兩種可能：1.2算法步驟〔1〕f1>f2,如圖1所示，極小點必在[a1,b]內，消去區(qū)間[a,a1),令a=a1，產生新區(qū)間[a,b],到此區(qū)間縮短了一次。值得注意的是新區(qū)間的a1點與原區(qū)間的a2點重合，可令a1=a2,這樣可少找一個新點和節(jié)省一次函數值計算?！?〕f1<=f2,極小點必在[a,a2]內，消去區(qū)間(a2,b],令b=a2,產生新區(qū)間[a,b],到此區(qū)間縮短了一次。同樣新區(qū)間a2點與原區(qū)間的a1點重合，可令a2=a1,f2<=f1?！?〕當縮短的新區(qū)間長度小于等于某一精度ε，即b-a<=ε時，取a*=(a1+a2)/2。1.3黃金分割法算法框圖給定a,b,e給定a,b,ea+0.382(b-a)a1a+0.382(b-a)a1,f(a1)f1a+0.618(b-a)a2,f(a2)f2NNYf1>f2?Yf1>f2?a=a1a=a1,a1=a2,f1=f2a2=a+0.618(b-a)f2=f(a2)b=a2,a2=a1,f2=f1a1=a+0.382(b-a)f1=f(a1)完畢輸出：a*=(a+b)/2NYb–a<e?2.拋物線插值法完畢輸出：a*=(a+b)/2NYb–a<e?2.1算法原理:拋物線也叫二次插值法，其理論依據為二次多項式可以在最優(yōu)點附近較好的逼近函數的形狀，做法是在函數f(x)的最優(yōu)點附近取三個構造點，然后用這三個點構造一條拋物線，把這條拋物線的極值點作為函數f(x)的極值點的近似。每次構造一條拋物線后，拋物線的極值點就可作為一個新的構造點，新的構造點與原來的三個構造點經過某種算法，得到下一步拋物線逼近的三個構造點，這就是拋物線法的算法過程。2.2算法步驟： 用拋物線求無約束問題minf(x),x∈R的算法步驟如下：①確定初始區(qū)間[a,b],選定初始插值內點t0∈(a,b)及精度e>0,令a0=a,b0=b,k=0;②求二次插值多項式的極小點：tk+1=.(2.10)③假設f(tk+1)≤f(tk),則當|tk+1-tk|≤e時，停頓迭代輸出tk+1;否則轉④;④假設f(tk+1)＞f(tk),則當|tk+1-tk|≤e時，停頓迭代輸出tk;否則轉⑤;⑤假設tk+1≤tk,令 ak+1=ak,bk+1=tk,tk+1=tk+1;置k=k+1轉②,否則令 ak+1=tk,bk+1=bk,tk+1=tk+1;置k=k+1轉②。⑥假設tk+1≤tk,令 ak+1=tk+1,bk+1=tk,tk+1=tk;置k=k+1轉②,否則令 ak+1=ak,bk+1=tk+1,tk+1=tk;置k=k+1轉②。初始插值內點可以取搜索區(qū)間[a,b]的中點。2.3拋物線插值法算法框圖開場開場確定t確定tk,ak,bk,要求f(ak)≥f(tk),f(bk)≥f(tk)按〔2.10〕計算tk+1按〔2.10〕計算tk+1t*=tk+1,f*=f(tk+1t*=tk+1,f*=f(tk+1)|tk+1-tk|＜eYN輸出t*,f*N輸出t*,f*Ntk+1＞Ntk+1＞tk?完畢完畢YYf(tf(tk+1)≤f(tk)?f(tk+1)≤f(tk)?f(tk+1)≤f(tk)?ak=tak=tk+1,f(ak)=f(tk+1)YYbk=tbk=tk,tk=tk+1,f(bk)=f(tk),f(tk)=f(tk+1)ak=tk,tk=tk+1,f(ak)=f(tk),f(tk)=f(tk+1)bk=tk+1,f(bk)=f(tk+1)3.算法的MATLAB實現3.1黃金分割法程序代碼在MATLAB中編程實現的黃金分割法函數為:golden。功能：用黃金分割法求解一維函數的極值。調用格式：tmin=golden(f,a,b,e)其中，f=目標函數； a=極值區(qū)間左端點；b=極值區(qū)間右端點； e=精度； tmin=目標取最小值時的自變量值；黃金分割法的MATLAB程序代碼如下：主程序：symstabe;a=input('搜索區(qū)間第一點\a=');b=input('搜索區(qū)間第二點\b=');e=input('搜索精度\ne=');disp('需求的最優(yōu)化函數f=f(t),tmin=golden(f,a,b,e)');m文件：functiontmin=golden(f,a,b,e)formatlong;k=0;a1=b-0.618*(b-a);a2=a+0.618*(b-a);whileb-a>ey1=subs(f,a1);y2=subs(f,a2);ify1>y2a=a1;a1=a2;y1=y2;a2=a+0.618*(b-a);elseb=a2;a2=a1;y2=y1;a1=b-0.618*(b-a);endk=k+1;endtmin=(a+b)/2;fmin=subs(f,tmin)fprintf('k=\n');disp(k);x=-3:0.001:5;y=x.*(x+2);plot(x,y,'k-',tmin,fmin,'bp');3.2實例驗證minf(t)=t*(t+2),初始單谷區(qū)間[a,b]=[-3,5],按精度e=0.001計算。代入到主程序代碼：搜索區(qū)間的第一點a=-3搜索區(qū)間的第二點b=5搜索精度e=0.001需求的最優(yōu)化函數f=f(t),tmin=golden(f,a,b,e)>>f=t*(t+2)f=t*(t+2)>>tmin=golden(f,a,b,e)k=19tmin=-1.000120312207862>>fmin=tmin*(tmin+2)fmin=-0.999999985524973目標函數的曲線圖如圖1-1所示:注：圖中☆所標識的點為黃金分割法所求的極小點。圖1-1函數f(t)=t(t+2)的曲線圖3.3誤差分析:1.f(t)=t*(t+2)的準確解如下：f=[1,2,0];p=polyder(f);t1=roots(p)f1=polyval(f,t1)t1=-1f1=-1極小點誤差：e1=(tmin-t1)/t1=(-1.000120312207862-(-1))/(-1)≈1e-4極小值誤差：e2=(f1-fmin)/f1=(-1+0.999999985524973)/(-1)≈1e-8結論：在精度e=0.001的情況下，通過黃金分割法經過19次迭代，求得：tmin=-1.000120312207862,與準確值的誤差e1=1e-4<e=0.001,滿足題目要求，黃金分割法是正確的，也可預見，隨著精度e的提高，tmin的值會越來越接近極小點，最終會等于準確值-1.由圖1-1也可以看出，圖中☆所標識的位置幾乎就在t=-1的位置上。3.4拋物線插值法程序代碼在MATLAB中編程實現的拋物線法函數為:minPWX。功能：用拋物線插值法求解一維函數的極值。調用格式：[x,minf]=minPWX(f,a,b,e)其中，f:目標函數；A:初始搜索區(qū)間左端點；b:初始搜索區(qū)間右端點； e:精度；x:目標取最小值時的自變量值；minf:目標函數的最小值。拋物線法的MATLAB程序代碼如下：function[x,minf]=minPWX(f,a,b,e)%目標函數：f;%初始搜索區(qū)間左端點：a;%初始搜索區(qū)間右端點：b;%精度：e;%目標函數取最小值時的自變量值:x;%目標函數的最小值:minfformatlong;ifnargin==3e=1.0e-3;endt0=(a+b)/2;k=0;tol=1;whiletol>efa=subs(f,findsym(f),a);%區(qū)間左端點函數值fb=subs(f,findsym(f),b);%區(qū)間右端點函數值ft0=subs(f,findsym(f),t0);%內插點函數值tu=fa*(b^2-t0^2)+fb*(t0^2-a^2)+ft0*(a^2-b^2);td=fa*(b-t0)+fb*(t0-a)+ft0*(a-b);t1=tu/2/td; %插值多項式的極小點ft1=subs(f,findsym(f),t1); %插值多項式的極小值tol=abs(t1-t0);ifft1<=ft0ift1<=t0b=t0;t0=t1;elsea=t0;t0=t1;endk=k+1;elseift1<=t0a=t1;elseb=t1;endk=k+1;endendx=t1;minf=subs(f,findsym(f),x);formatshort;3.5實例驗證用拋物線法求函數f(t)=t*(t+2)，初始單谷區(qū)間[a,b]=[-3,5],按精度e=0.001計算。解：在MATLAb命令窗口中輸入：>>symst;>>f=t*(t+2);>>