function nma_CG_TEST1() % driver function for nma_CG % % This function implements HW5 for Math 228A, UC Davis % Fall 2010. % It calls the CG() function written to solve conjugate % gradient using different preconditioners. % % After the call, the resulting table is plotted and the % tabled are formatted and displayed. % % by Nasser M. Abbasi % DEC 2010 % close all; myforce = @(X,Y) -exp( -(X-0.25).^2 - (Y-0.6).^2 ); %myforce = @(X,Y) exp(X).*sin(Y); %myforce = @(X,Y) 0.*X; methods = {'NONE' 'SSOR' 'MULTI-GRID' 'INCOMPLETE-CHOLESKY'}; methods = {'NONE'}; h = [1/16 1/32 1/64]; tol = 10^-6; max_iterations = 500; for k=1:length(methods) for i=1:length(h) n = (1/h(i) + 1); nx = n-2; %internal grid points [A,f] = nma_GENP2D(nx,myforce); t=cputime; tstart = tic; [u, nIter, result, cond_MA, cond_M, cond_A, eig_MA ,eig_M, eig_A] =... nma_CG( -A, 0, tol, max_iterations, methods{k},0,0.001 ); fprintf('CPU time for N=%d, method=%s is %f\n',nx^2, methods{k},cputime-t); telapsed = toc(tstart); fprintf('elapsed time for N=%d, method=%s is %f\n',nx^2, methods{k},telapsed); gen_tables(result,nIter,... sprintf('h =%s,n=%d\t eps=%f\tmethod=%s\n, k=%d, cond(1/M)*A=%f, cond A=%f, cond(M)=%f\n',... rats(h(i)),n,tol,methods{k},nIter,cond_MA,cond_A,cond_M)); % plot distribution of eigenvalues figure(); subplot(2,1,1); scatter(eig_A,zeros(length(eig_A),1),2,'r') t1=sprintf(... 'spectrum of A. $h=%s$, pre = %s, $max \\lambda=%f$, $min \\lambda=%f$',... rats(h(i)),methods{k},max(eig_A),min(eig_A)); t2=sprintf('condition number $=%f$',cond_A); title({t1;t2},'interpreter','latex'); xlabel('re( $\lambda$ )','interpreter','latex'); ylabel('im( $\lambda$ )','interpreter','latex'); subplot(2,1,2); scatter(eig_MA,zeros(length(eig_MA),1),2,'k') xlabel('re( $\lambda$ )','interpreter','latex'); ylabel('im( $\lambda$ )','interpreter','latex'); t1=sprintf(... 'spectrum of $M^{-1} A$, $max \\lambda=%f$, $min \\lambda=%f$',... max(eig_MA),min(eig_MA)); t2=sprintf('condition number $=%f$',cond_MA); title({t1;t2},'interpreter','latex'); % plot residual figure(); subplot(2,1,1); plot(log10(result(1:nIter,4)),'-o'); title(sprintf(... 'h=%s, pre-conditioner=%s, showing log10 of ||residual|| vs. iteration number',... rats(h(i)),methods{k})); xlabel('k, iteration number'); ylabel('log10(||r||)'); grid; subplot(2,1,2); plot(result(1:nIter,4),'-o'); title(sprintf(... 'h=%s, pre-conditioner=%s, showing ||residual|| vs. iteration number',... rats(h(i)),methods{k})); xlabel('k, iteration number'); ylabel('||r||'); grid; % plot error figure(); subplot(2,1,1); plot(log10(result(1:nIter,2)),'-o'); title(sprintf(... 'h=%s, pre-conditioner=%s, showing log10 of ||error|| vs. iteration number',... rats(h(i)),methods{k})); xlabel('k, iteration number'); ylabel('log10(||error||)'); grid; subplot(2,1,2); plot(result(1:nIter,2),'-o'); title(sprintf(... 'h=%s, pre-conditioner=%s, showing ||error|| vs. iteration number',... rats(h(i)),methods{k})); xlabel('k, iteration number'); ylabel('||error||'); grid; figure(); u=reshape(-u,nx,nx); mesh(u); title(sprintf('approximate solution found, h=%s, iterations=%d, method=%s',... rats(h(i)),nIter,methods{k})); end end end %-------------------------------- % function gen_tables(result,k,heading) titles={'k','|e|','ratio','|r|','ratio'}; wid = 16; fms = {'d','.7f','.7f','.7f','.7f'}; fileID=1; fprintf(fileID,heading); nma_format_matrix(titles,result(1:k,:),wid,fms,fileID,false); end