MATLAB的黄金分割法与抛物线插值法

...wd......wd......wd...MATLAB的黄金分割法与抛物线插值法摘要为了求解最优化模型的最优解，可使用基于MATLAB算法编程的黄金分割法与抛物线插值法，来实现求解的过程。黄金分割法是通过所选试点的函数值而逐步缩短单谷区间来搜索最优点，利用迭代进而得出结论。抛物线插值法亦称二次插值法，是一种多项式插值法，逐次以拟合的二次曲线的极小点，逼近原寻求函数极小点的一种方法。通过将MATLAB与最优化问题相结合，不仅可以加深对黄金分割法、抛物线插值法的根本理解和算法框图及其步骤的全面理解，也有利于帮助我们掌握MATLAB的使用方法。关键词：MATLAB，黄金分割法，抛物线插值法，最优解，迭代目录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引言为了对最优化问题进展求解，为了更好的掌握MATLAB的运用，本文主要介绍一维最优化问题的一维搜索法黄金分割法与抛物线插值法，利用MATLAB算法将其实现。黄金分割法也叫0.618法，它是一种基于区间收缩的极小点搜索算法，当用进退法确定搜索区间后，我们只知道极小点包含于搜索区间内，但是具体是哪个点，无法得知。黄金分割法的思想很直接，既然极小点包含于搜索区间内，那么可以不断的缩小搜索区间，就可以使搜索区间的端点逼近到极小点。抛物线插值法(parabolicinterpolationmet-hod)亦称二次插值法一种多项式插值法.逐次以拟合的二次曲线的极小点，逼近原寻求函数极小点的一种方法,能求得精度较高的解，但速度会比照慢。时至今日，经过MathWorks公司的不断完善，MATLAB已经开展成为适合多学科、多种工作平台的功能强劲的大型软件。在国外，MATLAB已经经受了多年考验。在欧美等高校，MATLAB已经成为线性代数、自动控制理论、数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的根本教学工具；成为攻读学位的大学生、硕士生、博士生必须掌握的根本技能。在设计研究单位和工业部门，MATLAB被广泛用于科学研究和解决各种具体问题。1.黄金分割法1.1算法原理黄金分割法的根本原理：黄金分割法又称0.618法，它是通过不断缩短搜索区间的长度来寻求一维函数的极小点。这种方法的原理是：在搜索区间[a,b]内按如下规则对称地取两点：计算它们的函数值f1=f(a1),f2=f(a2),比照它们的大小，结果有两种可能：1.2算法步骤〔1〕f1>f2,如图1所示，极小点必在[a1,b]内，消去区间[a,a1),令a=a1，产生新区间[a,b],到此区间缩短了一次。值得注意的是新区间的a1点与原区间的a2点重合，可令a1=a2,这样可少找一个新点和节省一次函数值计算。〔2〕f1<=f2,极小点必在[a,a2]内，消去区间(a2,b],令b=a2,产生新区间[a,b],到此区间缩短了一次。同样新区间a2点与原区间的a1点重合，可令a2=a1,f2<=f1。〔3〕当缩短的新区间长度小于等于某一精度ε，即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算法原理:抛物线也叫二次插值法，其理论依据为二次多项式可以在最优点附近较好的逼近函数的形状，做法是在函数f(x)的最优点附近取三个构造点，然后用这三个点构造一条抛物线，把这条抛物线的极值点作为函数f(x)的极值点的近似。每次构造一条抛物线后，抛物线的极值点就可作为一个新的构造点，新的构造点与原来的三个构造点经过某种算法，得到下一步抛物线逼近的三个构造点，这就是抛物线法的算法过程。2.2算法步骤： 用抛物线求无约束问题minf(x),x∈R的算法步骤如下：①确定初始区间[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转②。初始插值内点可以取搜索区间[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=极值区间左端点；b=极值区间右端点； e=精度； tmin=目标取最小值时的自变量值；黄金分割法的MATLAB程序代码如下：主程序：symstabe;a=input('搜索区间第一点\a=');b=input('搜索区间第二点\b=');e=input('搜索精度\ne=');disp('需求的最优化函数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),初始单谷区间[a,b]=[-3,5],按精度e=0.001计算。代入到主程序代码：搜索区间的第一点a=-3搜索区间的第二点b=5搜索精度e=0.001需求的最优化函数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:初始搜索区间左端点；b:初始搜索区间右端点； e:精度；x:目标取最小值时的自变量值；minf:目标函数的最小值。抛物线法的MATLAB程序代码如下：function[x,minf]=minPWX(f,a,b,e)%目标函数：f;%初始搜索区间左端点：a;%初始搜索区间右端点：b;%精度：e;%目标函数取最小值时的自变量值:x;%目标函数的最小值:minfformatlong;ifnargin==3e=1.0e-3;endt0=(a+b)/2;k=0;tol=1;whiletol>efa=subs(f,findsym(f),a);%区间左端点函数值fb=subs(f,findsym(f),b);%区间右端点函数值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)，初始单谷区间[a,b]=[-3,5],按精度e=0.001计算。解：在MATLAb命令窗口中输入：>>symst;>>f=t*(t+2);>>