牛顿法
迭代公式:(1)2()1()[()]()k k k k x x f x f x +-=-∇∇
Matlab 代码:
function [x1,k] =newton(x1,eps)
hs=inline('(x-1)^4+y^2'); 写入函数
ezcontour(hs,[-10 10 -10 10]); 建立坐标系
hold on; 显示图像
syms x y 定义变量
f=(x-1)^4+y^2; 定义函数
grad1=jacobian(f,[x,y]); 求f 的一阶梯度
grad2=jacobian(grad1,[x,y]); 求f 的二阶梯度
k=0; 迭代初始值
while 1 循环
grad1z=subs(subs(grad1,x,x1(1)),y,x1(2)); 给f 一阶梯度赋初值 grad2z=subs(subs(grad2,x,x1(1)),y,x1(2)); 给f 二阶梯度赋初值 x2=x1-inv(grad2z)*(grad1z)'; 核心迭代公式
if norm(x1-x2)<eps 判断收敛条件
break;
else
plot([x1(1),x2(1)],[x1(2),x2(2)],'-r*'); 画图
k=k+1; 迭代继续
x1=x2; 赋值
end
end
end
优点:在极小点附近收敛快
缺点:但是要计算目标函数的hesse 矩阵
最速下降法
1. :选取初始点xo ,给定误差
2. 计算一阶梯度。
若一阶梯度小于误差,停止迭代,输出
3. 取()()()k k p f x =∇
4. 10
t ()(), 1.min k k k k k k k k k k t f x t p f x tp x x t p k k +≥+=+=+=+进行一维搜索,求,使得令转第二步
例题:
求min (x-2)^4+(x-2*y)^2.初始值(0,3)误差为0.1
(1)编写一个目标函数,存为f.m
function z = f( x,y )
z=(x-2.0)^4+(x-2.0*y)^2;
end
(2)分别关于x 和y 求出一阶梯度,分别存为fx.m 和fy.m function z = fx( x,y )
z=2.0*x-4.0*y+4.0*(x-2.0)^3;
end
和
function z = fy( x,y )
z=8.0*y-4.0*x;
end
(3)下面是脚本文件,一维搜索用的是黄金分割法Tic 计算时间
eps=10^(-4);误差
err=10;
dt=0.01;
x0=1.0;初始值
y0=1.0;
mm=0;
while err>eps 黄金分割法
dfx=-fx(x0,y0);
dfy=-fy(x0,y0);
tl=0;tr=1;确定一维搜索的区间
h=3;
nn=0;
gerr=10;
geps=10^(-4);
while gerr>geps
tll=tl+0.382*abs(tr-tl);
trr=tl+0.618*abs(tr-tl);
if
f(x0+tll*h*dfx,y0+tll*h*dfy)>f(x0+trr*h*dfx,y0+trr*h*dfy) tl=tll;
else
tr=trr;
end
gerr=abs(tl-tr); 区间的长度之差
tt=0.5*(tl+tr);
nn=nn+1;步数增加
if nn>200 迭代终止条件
break
end
end
x0=x0+tt*h*dfx; 重新迭代
y0=y0+tt*h*dfy;
err=sqrt(fx(x0,y0)^2+fy(x0,y0)^2);
mm=mm+1;步数增加
if mm>700 迭代步数超过700,终止
break
end
end
res=[x0,y0];输出最后的x,y。
toc 计算运行时间
拟牛顿法(DFP 算法)
220'412010min ()4,(1,1),,1001f x x x x H ε-⎛⎫=+=== ⎪⎝⎭
取 这是一个脚本文件可以直接运行
syms x1 x2;定义变量
eps=0.00001;
x0=[1,1]';初始值
h0=[1,0;0,1];
f=x1^2+4*x2^2;待求函数
fx=diff(f,x1);对x求导
fy=diff(f,x2);对y求导
df=[fx,fy];f的一阶梯度
dfx0=[subs(fx,[x1,x2],x0),subs(fy,[x1,x2],x0)]';赋初值
d0=-dfx0;搜索方向
n=1;
while 1
syms t;
s0=x0+t*d0;引入变量t
ff=subs(f,[x1,x2],s0)给f赋值;
t=solve(diff(ff));求ff的极小点
xx1=x0+t*d0;更新初始值
dfx1=[subs(fx,[x1,x2],xx1'),subs(fy,[x1,x2],xx1')]';赋值
pp=sqrt(dfx1*dfx1');判断此时一阶梯度的值
if(pp<0.001)迭代终止条件
break
end
a1=xx1-x0;
r1=dfx1-dfx0;
h1=h0+(a1*a1')/(a1'*r1)-(h0*r1*r1'*h0)/(r1'*h0*r1);h0的更新d1=-h1*dfx1;搜索方向的更新
d0=d1;循环赋值
x0=xx1;循环赋值
h0=h1;二阶梯度的近似的更新
n=n+1;计算迭代步数
end
共轭梯度法
221212112(,)242f x x x x x x x =+--求二次函数的极小点。
(0)(1,1)T =取初始点x
这是一个脚本文件
clc;
clear all ;
syms x y t ; 定义变量
x0=[1,1];初始值
n=1;初始迭代
t=0;
f1=x^2+2*y^2-4*x-2*x*y;待求函数
dfx=diff(f1,x);求函数的对x 一阶梯度
dfy=diff(f1,y);函数对y 的一阶梯度
df=[dfx,dfy];函数一阶梯度以数组的形式
while 1
syms kk ;在循环里定义变量
g0=[subs(dfx,[x,y],x0),subs(dfy,[x,y],x0)];给一阶梯度赋值 s0=-g0;下降方向
m0=x0+kk*s0;引入变量kk
f11=f(m0(1),m0(2));带入原函数,得到关于kk的函数
kk=solve(diff(f11));求f11的极小点
m1=x0+kk*s0;更新迭代初始值
g1=[subs(dfx,[x,y],m1),subs(dfy,[x,y],m1)];给一阶梯度赋值 s1=-g1;
k=(s1*s1')/(g0*g0');
s2=s1+k*s0; 更新梯度
s0=s2;重新迭代
x0=m1;
tt=subs(df,[x,y],m1);
t=sqrt(tt*tt');一阶梯度值
n=n+1;
if (t<0.01)判断迭代终止条件
break
end
% if(n>20)
% break;
% end
end
最佳答案[4,2]。