工程数学上机实验四：共轭梯度法程序设计代码

function [k,x,val] = frcg(fun,gfun,x0,epsilon,N)

%共轭梯度法求解无约束问题

% fun,gfun分别为目标函数及其梯度，x0是初始点

% epsilon是容许误差，N是最大的迭代次数

if nargin<5, N=10000;end

if nargin<4, epsilon=1e-6;end

beta=0.6; sigma=0.4;

n=length(x0);

k=0;

while(k<N)

gk=feval(gfun,x0);%计算x0处的梯度方向

itern=k-(n+1)*floor(k/(n+1));

itern=itern+1;%每迭代n次后，要重新开始迭代

if(itern==1)

dk=-gk;%第一次的迭代方向取梯度的负方向

else

betak=(gk'*gk)/(g0'*g0);%gk'表示矩阵gk的共轭转置，计算β

dk=-gk+betak*d0;

gd=gk'*dk;

if(gd>=0.0)

dk=-gk;%当搜索方向不是梯度方向时，重新开始搜索

end

end

if(norm(gk)<epsilon)

break;

end

m=0;mk=0;

while(m<20)

if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk)

mk=m;

break;

end

m=m+1;

end

x=x0+beta^mk*dk;%计算下一个x

g0=gk;d0=dk;

x0=x;

k=k+1;%迭代次数加1

end

val=feval(fun,x);

end

 

fun.m文件

function f = fun(x)

f=(x(1)+10*x(2))^2+5*(x(3)-x(4))^2+(x(2)-2*x(3))^4+10*(x(1)-x(4))^4;

end

 

gfun.m文件

function gf = gfun(x)

gf=[2*(x(1)+10*x(2))+40*(x(1)-x(4))^3;20*(x(1)+10*x(2))+4*(x(2)-2*x(3))^3;10*(x(3)-x(4))-8*(x(2)-2*x(3))^3;-10*(x(3)-x(4))-40*(x(1)-x(4))^3];

end