# HG changeset patch # User Rik # Date 1373559954 25200 # Node ID b34202b2421201d05d6b9d40ce0472921fe32995 # Parent 1eb5e5f0ee133c03904fe87f67e9654d522bc3d4 fplot.m: Overhaul function for Matlab compatibility and performance (bug #38961). * scripts/plot/fplot.m: Add ability to specify n,tol,fmt in any order and simultaneously. Return data rather than plotting it if asked. Use additional test on progress of algorithm to decide whether to quit. Add %!demo and %!tests. diff -r 1eb5e5f0ee13 -r b34202b24212 scripts/plot/fplot.m --- a/scripts/plot/fplot.m Thu Jul 11 13:11:55 2013 +0100 +++ b/scripts/plot/fplot.m Thu Jul 11 09:25:54 2013 -0700 @@ -20,41 +20,51 @@ ## @deftypefn {Function File} {} fplot (@var{fn}, @var{limits}) ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{tol}) ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{n}) -## @deftypefnx {Function File} {} fplot (@dots{}, @var{fmt}) -## Plot a function @var{fn} within defined limits. +## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{fmt}) +## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{tol}, @var{n}, @var{fmt}) +## @deftypefnx {Function File} {[@var{x}, @var{y}] =} fplot (@dots{}) +## Plot a function @var{fn} within the range defined by @var{limits}. +## ## @var{fn} is a function handle, inline function, or string ## containing the name of the function to evaluate. -## The limits of the plot are given by @var{limits} of the form -## @code{[@var{xlo}, @var{xhi}]} or @code{[@var{xlo}, @var{xhi}, -## @var{ylo}, @var{yhi}]}. @var{tol} is the default tolerance to use for the -## plot, and if @var{tol} is an integer it is assumed that it defines the -## number points to use in the plot. The @var{fmt} argument is passed -## to the plot command. +## The limits of the plot are of the form @code{[@var{xlo}, @var{xhi}]} or +## @code{[@var{xlo}, @var{xhi}, @var{ylo}, @var{yhi}]}. +## The next three arguments are all optional and any number of them may be +## given in any order. +## @var{tol} is the relative tolerance to use for the plot and defaults +## to 2e-3 (.2%). +## @var{n} is the minimum number of of points to use. When @var{n} is +## specified, the maximum stepsize will be +## @code{@var{xhi} - @var{xlo} / @var{n}}. More than @var{n} points may still +## be used in order to meet the relative tolerance requirement. +## The @var{fmt} argument specifies the linestyle to be used by the plot +## command. +## +## With no output arguments the results are immediately plotted. With two +## output arguments the 2-D plot data is returned. The data can subsequently +## be plotted manually with @code{plot (@var{x}, @var{y}). +## +## Example: ## ## @example ## @group -## fplot ("cos", [0, 2*pi]) +## fplot (@cos, [0, 2*pi]) ## fplot ("[cos(x), sin(x)]", [0, 2*pi]) ## @end group ## @end example +## +## Note: @code{fplot} works best with continuous functions. Functions with +## discontinuities are unlikely to plot well. This restriction may be +## removed in the future. ## @seealso{plot} ## @end deftypefn ## Author: Paul Kienzle -function fplot (fn, limits, n = 0.002, fmt = "") - - if (nargin < 2 || nargin > 4) - print_usage (); - endif +function [X, Y] = fplot (fn, limits, varargin) - if (iscomplex (limits) || (numel (limits) != 2 && numel (limits) != 4)) - error ("fplot: LIMITS must be a real vector with 2 or 4 elements"); - endif - - if (ischar (n)) - fmt = n; - n = 0.002; + if (nargin < 2 || nargin > 5) + print_usage (); endif if (strcmp (typeinfo (fn), "inline function")) @@ -71,52 +81,89 @@ error ("fplot: FN must be a function handle, inline function, or string"); endif - if (n != fix (n)) - tol = n; - x0 = linspace (limits(1), limits(2), 5)'; + if (iscomplex (limits) || (numel (limits) != 2 && numel (limits) != 4)) + error ("fplot: LIMITS must be a real vector with 2 or 4 elements"); + endif + + n = 5; + tol = 2e-3; + fmt = ""; + for i = 1:numel (varargin) + arg = varargin{i}; + if (ischar (arg)) + fmt = arg; + elseif (isnumeric (arg) && isscalar (arg) && arg > 0) + if (arg == fix (arg)) + n = arg; + else + tol = arg; + endif + else + error ("fplot: bad input in position %d", i+2); + endif + endfor + + if (n != 5) + ## n was specified + x0 = linspace (limits(1), limits(2), n/2 + 1)'; y0 = feval (fn, x0); - err0 = Inf; - n = 8; x = linspace (limits(1), limits(2), n)'; y = feval (fn, x); - - if (! size_equal (x0, y0)) - ## FN is a constant value function - y0 = repmat (y0, size (x0)); - y = repmat (y, size (x)); - endif - - while (n < 2 .^ 20) - y00 = interp1 (x0, y0, x, "linear"); - err = 0.5 * max (abs ((y00 - y) ./ (y00 + y))(:)); - if (err == err0 || err < tol) - break; - endif - x0 = x; - y0 = y; - err0 = err; - n = 2 * (n - 1) + 1; - x = linspace (limits(1), limits(2), n)'; - y = feval (fn, x); - endwhile else + x0 = linspace (limits(1), limits(2), 5)'; + y0 = feval (fn, x0); + n = 8; x = linspace (limits(1), limits(2), n)'; y = feval (fn, x); endif - plot (x, y, fmt); - - if (length (limits) > 2) - axis (limits); + if (rows (x0) != rows (y0)) + ## FN is a constant value function + y0 = repmat (y0, size (x0)); + y = repmat (y, size (x)); endif - if (isvector (y)) - legend (nam); + err0 = Inf; + + ## FIXME: This algorithm should really use adaptive scaling as the + ## the numerical quadrature algorithms do so that extra points are + ## used where they are needed and not spread evenly over the entire + ## x-range. Try any function with a discontinuity such as + ## fplot (@tan, [-2, 2]) or fplot ("1./x", [-3, 2]) to see the + ## problems with the current solution. + + while (n < 2^18) # Something is wrong if we need more than 250K points + yi = interp1 (x0, y0, x, "linear"); + ## relative error calculation using average of [yi,y] as reference + ## since neither estimate is known a priori to be better than the other. + err = 0.5 * max (abs ((yi - y) ./ (yi + y))(:)); + if (err < tol || abs (err - err0) < tol/2) + ## Either relative tolerance has been met OR + ## algorithm has stopped making any reasonable progress per iteration. + break; + endif + x0 = x; + y0 = y; + err0 = err; + n = 2 * (n - 1) + 1; + x = linspace (limits(1), limits(2), n)'; + y = feval (fn, x); + endwhile + + if (nargout == 2) + X = x; + Y = y; else - for i = 1:columns (y) - nams{i} = sprintf ("%s(:,%i)", nam, i); - endfor - legend (nams{:}); + plot (x, y, fmt); + axis (limits); + if (isvector (y)) + legend (nam); + else + for i = 1:columns (y) + nams{i} = sprintf ("%s(:,%i)", nam, i); + endfor + legend (nams{:}); + endif endif endfunction @@ -124,17 +171,30 @@ %!demo %! clf; -%! fplot ('cos', [0, 2*pi]); +%! fplot (@cos, [0, 2*pi]); %!demo %! clf; %! fplot ('[cos(x), sin(x)]', [0, 2*pi]); +%!demo +%! clf; +%! ## sinc function +%! fh = @(x) sin (x*pi) ./ (x*pi); +%! fplot (fh, [-5, 5]); + +%!test +%! [x, y] = fplot ("[cos(x), sin(x)]", [0, 2*pi]); +%! assert (columns (y) == 2); +%! assert (rows (x) == rows (y)); +%! assert (y, [cos(x), sin(x)], -2e-3); + %% Test input validation %!error fplot (1) -%!error fplot (1,2,3,4,5) +%!error fplot (1,2,3,4,5,6) +%!error fplot (1, [0 1]) %!error fplot (@cos, [i, 2*i]) %!error fplot (@cos, [1]) %!error fplot (@cos, [1 2 3]) -%!error fplot (1, [0 1]) +%!error fplot (@cos,[-1,1], {1})