最近的暑期学校在上算法课,其中有一节课后作业便是用二次插值法求单变量函数在单峰区间上近似极小点,我需要自己编程设计算法程序。
算法书是马昌凤的**《最优化方法及其Matlab程序设计》**,本次编程也是依照该书的算法步骤,但是有改进。
- 算法简介
在求解一元函数f的极小点时,在搜索区间中用低次(通常不超过三次)插值多项式P来近似目标函数,后求该多项式的极小点(比较容易计算),并以此作为目标函数f的近似极小点。如果其近似的程度尚未达到所要求的精度时,反复使用此法,逐次拟合,直到满足给定的精度时为止。
-
算法原理及步骤
2.1
已知目标函数 ϕ \phi ϕ,为一元二次函数。
2.2 步骤
步0:给定区间[a,b],初始h=(b-a)/2,初始三点为:
s 0 = a , s 1 = s 0 + h , s 2 = s 0 + 2 ∗ h s0=a, s1=s0+h, s2=s0+2*h s0=a,s1=s0+h,s2=s0+2∗h
寻找合适的步长h, 确定三点s0<s1<s2,对应的函数值有:
φ 1 < φ 0 , φ 1 < φ 2 \qquad \qquad \qquad \qquad φ1<φ0, \quad φ1<φ2 φ1<φ0,φ1<φ2
在转步1前重新计算 s 1 = ( s 2 + s 0 ) / 2 和 f ( s 1 ) s1=(s2+s0)/2和f(s1) s1=(s2+s0)/2和f(s1)。 -
Matlab代码
function [s,k] = eccz(f,a,b,delt)
max=20;maxj=20;
k=0;j=0; %k为算法迭代次数,j为寻找合适h过程中迭代次数
h=(b-a)/2;
s0=a;s1=a+h;s2=a+2*h;
f0=feval(f,s0);
f1=feval(f,s1);
f2=feval(f,s2);
while(j<maxj && abs(h)>delt) %确定h使得f1<f0且f1<f2
if (f0<=f1)
s2=s1; f2=f1; h=0.5*h;
s1=s0+h;
f1=feval(f,s1);
else
if (f2<=f1)
s1=s2; f1=f2; h=2*h;
s2=s0+2*h;
f2=feval(f,s2);
end
end
j=j+1;
end
while(k<max && abs(s2-s0)>delt)
h=(s2-s0)/2;
hh=(3*f0-4*f1+f2)*0.5*h/(f0-2*f1+f2); %%由结论结果直接计算新h
h=hh;
ss=s0+hh;
fs=feval(f,ss);
if(fs>=f1)
if(s1<ss)
s2=ss;
f2=fs;
s1=(s2+s0)/2; %这里一定要计算一次,不然就会
f1=feval(f,s1); %出现程序一直循环不结束或者出现NAN情况!
else
s0=ss;
f0=fs;
s1=(s2+s0)/2;%同上
f1=feval(f,s1);
end
else
if(s1>ss)
s2=s1;
f2=f1;
s1=ss;
f1=fs;
else
s0=s1;
s1=ss;
f0=f1;
f1=fs;
end
end
k=k+1;
end
s=ss;
fprintf('经过%.1f次迭代,f的近似极小值点是:x=%.3f',k,s);
end
- 结果检验
以书上一道例题对本程序进行检验:
与题目答案吻合,但是迭代次数要比书上多6次。可惜书上的代码我看不太懂,应该不是按着书上算法步骤写的。