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