function ex8problem1 % Poisson matrix % -------------- n = 10; I = eye(n^2); A = gallery('poisson',n); % Choose a starting value %x0 = [1; zeros(n^2-1,1)]; x0 = ones(n^2,1); % Prolate matrix % -------------- % n = 10; % I = eye(n); % A = gallery('prolate',n); % Choose a starting value %x0 = [1; zeros(n-1,1)]; % x0 = ones(n,1); % Calculate exact solution: Xexact = min(eig(A)); % Define function handles here f = @(x) x'*(A*x)/(x'*x); df = @(x) 2*(I - x*x'/(x'*x))*(A*x)/(x'*x); % Line search parameters alpha0 = 1; beta = 0.5; c1 = 1e-4; tol = 1e-8; maxiter = 10000; %Run Fletcher-Reeves [X1,fX1,dfX1] = ncg(f,df,x0,c1,alpha0,beta,tol,maxiter,'fr'); %Run Polak-Ribiere [X2,fX2,dfX2] = ncg(f,df,x0,c1,alpha0,beta,tol,maxiter,'pr+'); %Run steepest descent [X3,fX3,dfX3] = steepdesc(f,df,x0,c1,alpha0,beta,tol,maxiter); %Plot function value and gradient subplot(1,2,1) semilogy([0:1:numel(fX1)-1],diag(dfX1'*dfX1).^(1/2),'ro'); hold on title('FR Gradient') subplot(1,2,2); semilogy([0:1:numel(fX1)-1],abs(fX1-Xexact),'ro'); hold on title('FR Error') subplot(1,2,1) semilogy([0:1:numel(fX2)-1],diag(dfX2'*dfX2).^(1/2),'ob'); title('PR+ Gradient') subplot(1,2,2); semilogy([0:1:numel(fX2)-1],abs(fX2-Xexact),'ob'); title('PR+ Error') subplot(1,2,1) semilogy([0:1:numel(fX3)-1],diag(dfX3'*dfX3).^(1/2),'ok'); title('SD') subplot(1,2,2); semilogy([0:1:numel(fX3)-1],abs(fX3-Xexact),'ok'); title('SD Error') end function [X,fX,dfX] = ncg(f,df,x0,c1,alpha0,beta,tol,maxiter,opt) X(:,1) = x0; fX(:,1) = f(x0); dfX(:,1) = df(x0); k = 1; xk = x0; gk = df(xk); pk = -gk; while norm(gk) > tol && k<=maxiter % -------------------------------- % ADD YOUR NONLINEAR CG CODE HERE % -------------------------------- end end function [X,fX,dfX] = steepdesc(f,df,x0,c1,alpha0,beta,tol,maxiter) X(:,1) = x0; fX(:,1) = f(x0); dfX(:,1) = df(x0); k = 1; xk = x0; gk = dfX(:,1); while norm(gk) > tol && k<=maxiter %search directions pk = -gk; % start backtracking alpha = alpha0; while f(xk + alpha*pk) > f(xk) + c1*alpha*gk'*pk alpha = alpha*beta; end %Perform step xk = xk + alpha*pk; gk = df(xk); X(:,k+1) = xk; fX(:,k+1) = f(xk); dfX(:,k+1) = gk; k = k+1; end end