function plotfuns(funname,a,b,tol,details) % Help for plotfuns.m % copyright 1998, D. Estep % % PURPOSE % % This script plots user input functions on [a,b] on one plot at points determined % using an adaptive strategy that insures the change in each function's % values from one sample point to the next is smaller than a tolerance. % The program starts with a nearly uniform mesh of sample points and then % creates an adapted mesh by successive refinement. A different mesh is % computed for each function. % % FUNCTION STATEMENT % % plotfuns(funname,a,b,tol,details) % % REQUIRED VARIABLES % % funname: contains the name of the functions to be plotted, the functions % can be built-in matlab functions and/or defined as m files % % a,b: numbers defining the interval [a,b] on which the function is plotted % % OPTIONAL VARIABLES % % tol: the tolerance for the maximum possible change in the function % values at neighboring sample points. % details: indicates whether the program should give information on % the sample points. % details = 0 means just show plot % details = 1 means show plot, changes in function, and mesh sizes % % USAGE % % First define a variable that contains the names of the functions to be plotted % using the matlab char statement. For instructions on using char, see help char. % The variable is defined using a matlab statement of the form % >> variable name = char('function1(x)','function2(x)', ...) % where function1, function1, ... are the names of the functions, which must % be functions of x. % % For example, we would define % >> funname = char('sin(x)','x*exp(x)') % if we wanted to plot sin(x) and x*exp(x) on the same plot. % % A function can also be defined in a m file. To do this, create a m % file named functionname.m containing the function to be plotted % with the form: % % function [y] = functionname(x) % y = give the formula in terms of input x; % % For example, we create a file called quadratic.m containing the lines % function [y] = quadratic(x) % y = x^2; % % Now when defining the variable containing the function name, we would use % 'functionname(x)' as one of the entries. % % For the example above, we would define % >> funname = char('sin(x),'quadratic(x)') % to get a plot of the matlab function sin(x) and the user defined % function quadratic(x) on the same plot. % % Finally, define the values of a, b, tol, details, then call % >> plotfuns(variable containing function names, a, b, tol, details) % % In the example above, we would type % >> plotfuns(funname,-4,4,.01,0) % to get a plot on the interval [-4,4] using a mesh of sample points % in which the change in the function values is smaller than .01 % % The last two arguments to plotfun can be left off, for example % >> plotfuns(funname,-4,4) % % % DETAILS % % The program uses two meshes: a current mesh and a new mesh. Going from % left to right, the difference in function values at successive nodes in % the current mesh is checked. If this difference is smaller than tol, the % new mesh points are defined to be the old mesh points. If this difference % is larger than tol, then the interval between the nodes is divided into % uniform pieces and the new mesh is defined by using the nodes of the % pieces. The number of pieces is computed by assuming the function is % linear between the nodes of the current mesh. After every interval in % the current mesh is checked, the current mesh is reset to be the just % computed new mesh. If no refinements were needed, the iteration stops % and if there was a refinement in some interval, there is another loop. %initialize % define the maximum number of adaptive iterations for any one function maxrefineiter=5; % check the number of arguments and set defaults for the missing values noargs = nargin; if noargs == 3 tol = .1 *(b-a); details = 0; elseif noargs == 4 details = 0; elseif noargs < 3 & noargs>5 disp('wrong number of inputs!') break end % nf is the number of functions to be plotted [nf lf] = size(funname); % loop over the functions for jf = 1: nf fprintf('computing plot points for function number %d\n',jf) % define initial number of sample points and uniform sample interval n(jf) = 11; delta = (b-a)/(n(jf)-1); % define the function to be plotted, the function name is in the jf row of funname onefunname = funname(jf,1:lf); f = inline(onefunname,'x'); % define the sample points. The first and last points align with a and b. x(1,jf) = a; y(1,jf) = f(x(1,jf)); x(n,jf) = b; y(n,jf)=f(x(n(jf),jf)); for i = 2: n(jf)-1 x(i,jf) = a+ delta*(i-1); % vary the internal sample points by 10% randomly to avoid % problems with sampling at symmetric points x(i,jf) = x(i,jf) + (rand-.5)*.2*delta; y(i,jf)=f(x(i,jf)); end % set refinement level refine = -1; % loopflag=1 means do a refinement search, the first time we do one of course loopflag=1; while loopflag>0 & refine <= maxrefineiter % xnew, ynew are the new mesh points and function values % nnew is the number of mesh points in the new mesh. These % are temporary variables and are not saved. % the first nodes in two meshes always agree xnew(1) = x(1,jf); ynew(1) = y(1,jf); nnew=1; % default: no more refinement iterations loopflag = -1; refine = refine+1; for i = 2:n(jf) % check if function difference is small enough if (abs(y(i,jf)-y(i-1,jf))< tol) % if so, then set new mesh point = old mesh point and go to next nnew=nnew+1; xnew(nnew) = x(i,jf); ynew(nnew) = y(i,jf); else % indicate that a refinement is needed in this iteration, so check % for one more iteration to be sure that the differences are small loopflag = 1; % assume that the function is linear between the two mesh points and % compute the necessary number of divisions to get the desired change nmid=ceil(abs(y(i,jf)-y(i-1,jf))/tol)+1; delta = abs(x(i,jf)-x(i-1,jf))/nmid; % now define the new mesh points up to the next to last one for k = 1:nmid-1; xnew(nnew+k) = x(i-1,jf) + delta*k; ynew(nnew+k) = f(xnew(nnew+k)); end % treat the last point to be sure that it lines up with old mesh xnew(nnew+nmid) = x(i,jf); ynew(nnew+nmid) = y(i,jf); nnew = nnew + nmid; end end % prepare for a new loop, set the old values to be the new values % this works because xnew and ynew are always longer or have the % same length as x and y. for i = 1: nnew x(i,jf) = xnew(i); y(i,jf) = ynew(i); end n(jf)=nnew; clear xnew; clear ynew; clear nnew; end fprintf('plotting took %d refinement steps\n', refine) if refine > maxrefineiter fprintf('The program exceeded %d refinement iterations plotting function number %d \n',maxrefineiter,jf) return end end % Because the function use different numbers of plotting points, some of % the entries at the ends of the columns of x and y may be undefined. % Fill in these entries with the last value that is defined. M=norm(n,inf); for jf = 1: nf for i = n(jf):M x(i,jf)=x(n(jf),jf); y(i,jf)=y(n(jf),jf); end end if(details==0) plot(x,y) %graphof=strcat('Graph of ',funname); %title('Plots of Functions') xlabel('x') else figure subplot(1,3,1) hold on plot(x,y) %plot(x,0,'r.') %graphof=['Graph of ',funname]; title('Plots') xlabel('x') hold off for jf = 1: nf for i = 1: M - 1 dy(i,jf) = abs(y(i+1,jf)-y(i,jf)); dx(i,jf) = abs(x(i+1,jf)-x(i,jf)); xx(i,jf) = (x(i+1,jf)+x(i,jf))/2; end end subplot(1,3,2) plot(xx,dy) title('Changes in Values') xlabel('x') subplot(1,3,3) semilogy(xx,dx) title('Interval Sizes') xlabel('x') end