以下程序能够适应一维和多维的目标函数。
powell.m
function [x,fval]=powell(f,x0,e,n)
% Powell 函数f 初始点x0 精度e 未知数个数n
X=zeros(n,n+2); % X矩阵第一列为初始点,第n+1列为Sn方向搜索后的点,最后一列为反射点
X(:,1)=x0; % 将初始点赋值给矩阵X
S=eye(n); % n阶单位矩阵的n列向量分别代表n元方程的n个搜索方向
delta=zeros(1,n); % 预分配内存,用于存放各点的函数下降量
x=x0';
fval=f(x0);
k=0;
while(k<20) % 最大迭代次数暂定为20
k=k+1;
for i=1:n
syms a;
a=solve(diff(subs(f(X(:,i)+a*S(:,i))),a)==0); % 求得当前搜索方向下最短步长
X(:,i+1)=X(:,i)+a*S(:,i); % 迭代后的新点坐标
delta(i)=f(X(:,i))-f(X(:,i+1)); % 函数值下降量
end
% 判断是否符合精度,退出循环
x=X(:,1)';
fval=f(X(:,1));
if(norm(X(:,n+1)-X(:,1))<=e)
break
end
S(:,n+1)=X(:,n+1)-X(:,1); % 生成新的搜索