japprox.m

% Example: japprox
% ~~~~~~~~~~~~~~~~
% This program approximates the value for the
% polar moment of inertia, J, for a circular
% cross section.  The circle is approximated
% using an n-sided polygon.
%
% Data is defined in the declaration statements
% below, where:
%
% radius  - radius of circle with center
%           located at (0,0)
% epsilon - allowable error tolerance 
%           (epsilon*100 gives the maximum
%           percent error allowed)
% maxarcs - maximum number of arcs to attempt
%           to subdivide circle into trying
%           to satisfy the error tolerance
%
% NOTES: Requires-
%           37 arcs for max of 1.00% error
%          115 arcs for max of 0.10% error
%          363 arcs for max of 0.01% error
%
% User m functions required:
%    circle, genprint
%----------------------------------------------

clear;
%...Input definitions
radius=5; epsilon=1e-2; maxarcs=400;

Japprox=10e20;  % force thru loop first time
No_arcs=4;      % square is first approximation

J=pi*radius^4/2; % true value of J

i=0; 
while abs((J-Japprox)/J) > epsilon
  fprintf('.');
  i=i+1; 
  if i > maxarcs
    fprintf('\n\nMaximum subdivision exceeded');
    fprintf('\nArcs attempted: %g\n\n',xp(i-1));
    error('Program exit');
  end
  No_arcs=No_arcs+1;
  [x,y]=circle(No_arcs,0,0,radius);
  [A,Qx,Qy,Ix,Iy,Ixy,xbar,ybar]=prop(x,y);
  Japprox=Ix+Iy; Jcalc(i)=Japprox;
  xp(i)=No_arcs;
end
last=length(xp);
PerErr=abs((J-Jcalc(last))/J)*100;

fprintf('\n\nApproximation of J');
fprintf(    ' for a Circle');
fprintf(  '\n------------------');
fprintf(    '-------------');
fprintf('\n\nEpsilon:        %g',epsilon);
fprintf(  ' or %g%%',epsilon*100);
fprintf(  '\nArcs required:  %g',xp(last));
fprintf(  '\nPercent error:  %g%%',PerErr);
fprintf(  '\nTrue J:         %g',J);
fprintf(  '\nApproximate J:  %g',Jcalc(last));
fprintf('\n');

%...Plot
tit=['Approximation of J for Circle',...
     ' (Maximum Percent Error=', ...
     num2str(epsilon*100),'%)'];
clf; 
  plot([min(xp) max(xp)],[J J],'-');
  hold on; plot(xp,Jcalc,'--');
  title(tit); xlabel('Number of Arcs Used');
  ylabel('J (L^4)');
  tmp=legend(' J',' Approximate J','Location','Best');
  axes(tmp); hold off; drawnow;
% genprint('geomerr');