matlab使用sym类型，subs和diff函数，进行模式搜索法和改进的powell法求解函数

如果想以分数输出结果，可以用   format rat

f函数使用以前的稍加改造

function y=f(x)
if(length(x)==1)
    global xk;
    global pk;
    x=xk+x*pk;
end
y=(-x(1)+x(2)+x(3))^2+(x(1)-x(2)+x(3))^2+(x(1)+x(2)-x(3))^2;

1.模式搜索法求解

初始点为（1，1）T

clc
clear
global xk;
global pk;
syms x1 x2 a;
y=x1^2+2*x2^2-4*x1-2*x1*x2;
fprintf('用模式搜索法求解问题:\n');
fprintf('min f(x)=x1^2+2*x2^2-4*x1-2*x1*x2\n');
fprintf('初始点为  \n');
%step 0
x=[1;1]
e=0.001;
n=2
p=eye(n);    %二阶单位矩阵
k=1;
while 1
    fprintf('=======================================\n');
    fprintf('第%d次:\n',k);
    %step 1
    xk=x;
    i=1;
    while 1
        fprintf('i=%d:-------------------------------\n',i);
        pk=p(:,i)
        y1=subs(subs(y,x1,xk(1)+a*pk(1)),x2,xk(2)+a*pk(2));
        fprintf('f(x+a*pk)  = ');
        disp(y1);
        y1=diff(y1);
        fprintf('f’(x+a*pk) = ');
        disp(y1);
        a1=solve(y1);
        a1=eval(a1);
        fprintf('a = ');
        disp(a1);
        xk=xk+a1*pk;
        fprintf('xk = \n');
        disp(xk);
        fprintf('-------------------------------\n');
        %step 3
        if(i<n)
            i=i+1;
        else
            break;
        end
    end
    %step 4
    %范数用的是平方和开根号
    xk_x=sum((xk-x).^2)
    xk_x=sqrt(xk_x)
    if xk_x<=e
        fprintf('x* = \n');
        disp(x);
        break;
    else
        pk=xk-x
        y1=subs(subs(y,x1,xk(1)+a*pk(1)),x2,xk(2)+a*pk(2));
        fprintf('f(x+a*pk)\t= ');
        disp(y1);
        y1=diff(y1);
        fprintf('f’(x+a*pk)\t= ');
        disp(y1);
        a1=solve(y1);
        a1=eval(a1);
        fprintf('a = ');
        disp(a1);
        xk=xk+a1*pk;
        fprintf('xk = \n');
        disp(xk);
    end
    x=xk;
    k=k+1;
end

2.用改进的powell算法求解

初始点为（1/2,1,1/2）

clc
clear
global xk;
global pk;
syms x1 x2 x3 a;
y=(-x1+x2+x3)^2+(x1-x2+x3)^2+(x1+x2-x3)^2;
fprintf('用改进的Powell法求解问题:\n');
fprintf('min f(x)=(-x1+x2+x3)^2+(x1-x2+x3)^2+(x1+x2-x3)^2\n');
n=3
fprintf('初始点为  \n');
%step 1
x=[1/2,0,0,0;1,0,0,0;1/2,0,0,0];
disp(x(:,1))
fx=[0,0,0,0];
fx(1)=f(x(:,1));
e=0.001;
p=eye(3);    %三阶单位矩阵
c=1;
k=1;
while 1
    fprintf('=======================================\n');
    fprintf('第%d次:\n',c);
    %step 2
    while 1
        fprintf('k=%d:-------------------------------\n',k);
        pk=p(:,k);
        xk=x(:,k)
        y1=subs(subs(subs(y,x1,xk(1)+a*pk(1)),x2,xk(2)+a*pk(2)),x3,xk(3)+a*pk(3));
        fprintf('f(x+a*pk) = ');
        disp(y1);
        y1=diff(y1);
        fprintf('f’(x+a*pk) = ');
        disp(y1);
        a1=solve(y1);
        a1=eval(a1);
        fprintf('a = ');
        disp(a1);
        x(:,k+1)=x(:,k)+a1*pk;
        fprintf('xk = \n');
        disp(x(:,k+1));
        fprintf('-------------------------------\n');
        %step 3
        if(k<n)
            k=k+1;
        else
            break;
        end
    end
    %step 4
    %范数用的是平方和开根号
    xn_x0=sqrt(sum((x(:,n+1)-x(:,1)).^2))
    if sqrt(sum((xn_x0.^2)))<=e
        fprintf('x* = \n');
        disp(x(:,n+1));
        break;
    end
    %setp 5
    k=1;
    m=1;
    fx(k+1)=f(x(:,k+1));
    d=fx(k)-fx(k+1);
    k=k+1;
    while 1
        fx(k+1)=f(x(:,k+1));
        d_=fx(k)-fx(k+1);
        if(d_>d)
            d=d_;
            m=k;
        k=k+1
        end
        if(k==n)
            break;
        end
    end
    f_=f(2*x(:,n+1)-x(:,1))
    %step 6
    if (f_>=fx(1))
        fx(1)=fx(n+1);
        x(:,1)=x(:,n+1);
        k=1;
    else
        if (fx(1)-2*fx(n+1)+f_)*(fx(1)-fx(n+1)-d)^2>=0.5*(fx(1)-f_)^2*d
            fx(1)=fx(n+1);
            x(:,1)=x(:,n+1);
            k=1;
        else
            while 1
                p(:,m)=p(:,m+1);
                m=m+1;
                if(m==n)
                    break;
                end
            end
            k=1;
            pk=p(:,n);
            xk=x(:,n+1)
            p(n)=(x(:,n+1)-x(:,1))/xn_x0;
            y1=subs(subs(subs(y,x1,xk(1)+a*pk(1)),x2,xk(2)+a*pk(2)),x3,xk(3)+a*pk(3));
            fprintf('f(x+a*pk) = ');
            disp(y1);
            y1=diff(y1);
            fprintf('f’(x+a*pk) = ');
            disp(y1);
            a1=solve(y1);
            a1=eval(a1);
            fprintf('a = ');
            disp(a1);
            x(:,1)=x(:,n)+a1*pk;
            fprintf('xk = \n');
            disp(x(:,1));
        end
    end
    c=c+1;
end