Mercurial > hg > octave-nkf
changeset 16950:b34202b24212
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.
author | Rik <rik@octave.org> |
---|---|
date | Thu, 11 Jul 2013 09:25:54 -0700 |
parents | 1eb5e5f0ee13 |
children | a7c9be4a2c0f |
files | scripts/plot/fplot.m |
diffstat | 1 files changed, 119 insertions(+), 59 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/plot/fplot.m +++ b/scripts/plot/fplot.m @@ -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 <pkienzle@users.sf.net> -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 <FN must be a function handle> fplot (1, [0 1]) %!error <LIMITS must be a real vector> fplot (@cos, [i, 2*i]) %!error <LIMITS must be a real vector with 2 or 4> fplot (@cos, [1]) %!error <LIMITS must be a real vector with 2 or 4> fplot (@cos, [1 2 3]) -%!error <FN must be a function handle> fplot (1, [0 1]) +%!error <bad input in position 3> fplot (@cos,[-1,1], {1})