最优化牛顿法最速下降法共轭梯度法matlab代码

牛顿法

迭代公式：(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)

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;循环赋值