共轭梯度法的MATLAB实现程序 共轭梯度法的 实现程序 使用MATLAB编程实现共轭梯度法 %conjugate gradient methods %method:FR,PRP,HS,DY,CD,WYL,LS %精确线搜索,梯度终止准则 function [ m,k,d,a,X,g1,fv] = conjgradme( G,b,c,X,e,method) if nargin<6 error('输入参数必须为6'); end n=length(G); if n==2 format long e %rat syms x1 x2 f=1/2*[x1,x2]*G*[x1;x2]+b'*[x1;x2]+c; g=[diff(f,x1);diff(f,x2)]; g1=subs(subs(g,x1,X(1,1)),x2,X(2,1)); d=-g1; a=-(d'*g1)/(d'*G*d);% a=-((X(:,1)'*G*d+b'*d)/(d'*G*d)); a=g1(:,1)'*g1(:,1)/(d(:,1)'*G*d(:,1)); X(:,2)=X(:,1)+a*d; g1=[g1 subs(subs(g,x1,X(1,2)),x2,X(2,2))]; m1=norm(g1(:,1)); m=norm(g1(:,2)); i=2; k=zeros(1); switch method case 'FR' while m>=e k(i-1)=(m/m1)^2; d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i)); %a1(i)=-((X(:,i)'*G*d(:,i)+b'*d(:,i))/(d(:,i)'*G*d(:,i)));a(i)=g1(:,i)'*g1(:,i)/(d(:,i)'*G*d(:,i)); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))]; m1=m; m=norm(g1(:,i+1)); i=i+1; end 

case 'PRP' while m>=e k(i-1)=g1(:,i)'*(g1(:,i)-g1(:,i-1))/(norm(g1(:,i-1)))^2; d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i)); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))]; m=norm(g1(:,i+1)); i=i+1; end case 'HS' while m>=e k(i-1)=g1(:,i)'*(g1(:,i)-g1(:,i-1))/(d(:,i-1)'*(g1(:,i)-g1(:,i-1))); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i)); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))]; m=norm(g1(:,i+1)); i=i+1; end case 'DY' while m>=e k(i-1)=g1(:,i)'*g1(:,i)/(d(:,i-1)'*(g1(:,i)-g1(:,i-1))); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i)); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))]; m=norm(g1(:,i+1)); i=i+1; end case 'LS' while m>=e k(i-1)=g1(:,i)'*(g1(:,i)-g1(:,i-1))/(d(:,i-1)'*(-g1(:,i-1))); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i)); %a(i)=-((X(:,i)'*G*d(:,i)+b'*d(:,i))/(d(:,i)'*G*d(:,i))); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))]; 

m=norm(g1(:,i+1)); i=i+1; end case 'CD' while m>=e k(i-1)=g1(:,i)'*g1(:,i)/(d(:,i-1)'*(-g1(:,i-1))); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i)); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))]; m=norm(g1(:,i+1)); i=i+1; end case 'WYL' while m>=e k(i-1)=g1(:,i)'*(g1(:,i)-(m/m1)*g1(:,i-1))/(m1^2); d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1); a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i)); %a(i)=-((X(:,i)'*G*d(:,i)+b'*d(:,i))/(d(:,i)'*G*d(:,i))); X(:,i+1)=X(:,i)+a(i)*d(:,i); g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))]; m1=m; m=norm(g1(:,i+1)); i=i+1; end end fv=subs(subs(f,x1,X(1,i)),x2,X(2,i)); end l1=X(1,i);l2=X(2,i); w1=X(1,1);w2=X(2,1); v1=min(l1,w1)-abs(l1-w1)/10:abs(l1-w1)/10:max(l1,w1)+abs(l1-w1)/10; v2=min(l2,w2)-abs(l2-w2)/10:abs(l2-w2)/10:max(l2,w2)+abs(l2-w2)/10; [x,y]=meshgrid(v1,v2); s=size(x); z=zeros(size(x)); for i=1:s(1) for j=1:s(2) z(i,j)=1/2*[x(i,j),y(i,j)]*G*[x(i,j);y(i,j)]+b'*[x(i,j);y(i,j)]+c; 

end end [px,py] = gradient(z,.2,.2); contour(v1,v2,z), hold on, quiver(v1,v2,px,py) [C,h] = contour(x,y,z); set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2) x1=X(1,:);y1=X(2,:); plot(x1,y1,'r*:');