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