查看: 410|回复: 0|关注: 0

[未答复] 利用matlab实现基于阻尼牛顿法函数的极小值确定

[复制链接]

新手

5 麦片

财富积分


050


2

主题

3

帖子

0

最佳答案
发表于 2019-9-12 21:24:23 | 显示全部楼层 |阅读模式
% ******forwardBackward.m******%
function [minx,maxx] = forwardBackward (f,x0,h,eps)
%目标函数:f;
%初始点:x0;
% 初始步长:h;
% 精度:eps;
% 目标函数读取包含极值区间左端点:minx;
% 目标函数读取包含极值区间右端点:maxx;
% 处理输入三个参数时情况
if nargin == 3
  eps = 1.0e-6;
end
x1 = x0;
f1 = subs(f,findsym(f),x1);
x2 = x0 + h;
f2 = subs(f,findsym(f),x2);
if f1 > f2
  x1 = x2;
  x3 = x2;
  x2 = x1 + h;
  f1 = subs(f,findsym(f),x1);
  f2 = subs(f,findsym(f),x2);
  while f1 > f2
h = 2 * h;
x1 = x3;
x2 = x1 + h;
f1 = subs(f,findsym(f),x1);
f2 = subs(f,findsym(f),x2);
  end
  minx = x1;
  maxx = x2;
else if f1 == f2
minx = x1;
maxx = x2;
    else  
h = -h;
x2 = x1;
x3=x1;
x1 = x2 + h;
f1 = subs(f,findsym(f),x1);
f2 = subs(f,findsym(f),x2);
while f1 < f2
h = 2 * h;
x2 = x3;
x1 = x2 + h;
f1 = subs(f,findsym(f),x1);
f2 = subs(f,findsym(f),x2);
end
minx = x1;
maxx = x2;
end
end

%******quadraticInterpolation.m******%
function [ minx,minf ] = quadraticInterpolation( f,x1,x3,eps1,eps2 )
%目标函数:f;
%区间上界:a;
% 区间下界:b;
% 精度:eps;
% 目标函数极值点:minx;
% 处理输入三个参数时情况
if nargin ==3
    eps1 =1.0e-6;
    eps2 =1.0e-6;
end
x2=(x1+x3)/2;
i=1;
while i==1
    f1=feval(f,x1);
    f2=feval(f,x2);
    f3=feval(f,x3);
    k1=(x2-x3)*f1+(x2-x1)*f2+(x1-x2)*f3;
    k2=(x2^2-x3^2)*f1+(x3^2-x1^2)*f2+(x1^2-x2^2)*f3;
    if abs(k1)<eps1
        minx=x2;
        minf=feval(f,x2);
        break;
    else
        xp=k2/(2*k1);
        c1=xp-x1;
        c2=x3-xp;
        if c1*c2<=0
            minx=x2;
            minf=feval(f,x2);
            break;
        else
            if abs(xp-x2)<eps2
                if feval(f,xp)<f2
                    minx=xp;
                    minf=feval(f,xp);
                else
                    minx=x2;
                    minf=feval(f,x2);
                end
                break;
            else
                if xp>x2
                    if f2>feval(f,xp)
                        x1=x2;
                        x3=xp;
                        x2=0.5*(x1+x3);
                    else
                        x3=xp;
                    end
                end

                if f2==feval(f,xp)
                    x2=xp;
                    x3=x2;
                else
                    if f2==feval(f,xp)
                        x1=xp;
                        x3=x2;
                        x2=0.5*(x1+x3);
                    else
                        x1=xp;
                    end
                end
            end
        end
    end            
end

%******grampNewTonOptimizMethod.m******%
clc
clear
eps=0.001;%精度
x0=[5;0.5];%设置搜索起始点`
syms x1 x2 fx1 fx2 alpha;
f='x1^2+5*x2^2-2*x1*x2';
objectFun=inline(f,'x1','x2');
fx1=diff(f,x1);
fx2=diff(f,x2);
f2x1=diff(fx1,x1);
f2x2=diff(fx2,x2);
fx1x2=diff(fx1,x2);
G=[f2x1 fx1x2;fx1x2 f2x2];%获得海森矩阵
gx1=inline(fx1);
gx2=inline(fx2);
inverseG=inv(G);
gx1Value=feval(gx1,x0(1),x0(2));
gx2Value=feval(gx2,x0(1),x0(2));
deltf=[feval(gx1,x0(1),x0(2));feval(gx2,x0(1),x0(2))];%当前迭代点梯度
modelf=norm(deltf);
while modelf>eps
    x=x0-alpha*inverseG*deltf;
    falpha=feval(objectFun,x(1),x(2));
    falpha=inline(falpha);%目标函数转化为@函数
    [minx,maxx]=forwardBackward(falpha,1,0.1);
    t=quadraticInterpolation(falpha,minx,maxx);
    x=x0-t*inverseG*deltf;
    x0=x;
    deltf=[feval(gx1,x0(1),x0(2));feval(gx2,x0(1),x0(2))];
    modelf=norm(subs(deltf));
end
  objectPoint=subs(x0)
  fmin=subs(feval(objectFun,x0(1),x0(2)))

得出结果为
objectPoint =
   5
1/2
fmin =
85/4
和给定的答案不一样,不知道哪里有问题。

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表