Mercurial > hg > octave-nkf

comparison scripts/plot/fplot.m @ 16950:b34202b24212

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 5ec3f4aea91c
children a7c9be4a2c0f
comparison
equal deleted inserted replaced
16949:1eb5e5f0ee13 16950:b34202b24212
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {} fplot (@var{fn}, @var{limits}) 20 ## @deftypefn {Function File} {} fplot (@var{fn}, @var{limits})
21 ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{tol}) 21 ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{tol})
22 ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{n}) 22 ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{n})
23 ## @deftypefnx {Function File} {} fplot (@dots{}, @var{fmt}) 23 ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{fmt})
24 ## Plot a function @var{fn} within defined limits. 24 ## @deftypefnx {Function File} {} fplot (@var{fn}, @var{limits}, @var{tol}, @var{n}, @var{fmt})
25 ## @deftypefnx {Function File} {[@var{x}, @var{y}] =} fplot (@dots{})
26 ## Plot a function @var{fn} within the range defined by @var{limits}.
27 ##
25 ## @var{fn} is a function handle, inline function, or string 28 ## @var{fn} is a function handle, inline function, or string
26 ## containing the name of the function to evaluate. 29 ## containing the name of the function to evaluate.
27 ## The limits of the plot are given by @var{limits} of the form 30 ## The limits of the plot are of the form @code{[@var{xlo}, @var{xhi}]} or
28 ## @code{[@var{xlo}, @var{xhi}]} or @code{[@var{xlo}, @var{xhi}, 31 ## @code{[@var{xlo}, @var{xhi}, @var{ylo}, @var{yhi}]}.
29 ## @var{ylo}, @var{yhi}]}. @var{tol} is the default tolerance to use for the 32 ## The next three arguments are all optional and any number of them may be
30 ## plot, and if @var{tol} is an integer it is assumed that it defines the 33 ## given in any order.
31 ## number points to use in the plot. The @var{fmt} argument is passed 34 ## @var{tol} is the relative tolerance to use for the plot and defaults
32 ## to the plot command. 35 ## to 2e-3 (.2%).
36 ## @var{n} is the minimum number of of points to use. When @var{n} is
37 ## specified, the maximum stepsize will be
38 ## @code{@var{xhi} - @var{xlo} / @var{n}}. More than @var{n} points may still
39 ## be used in order to meet the relative tolerance requirement.
40 ## The @var{fmt} argument specifies the linestyle to be used by the plot
41 ## command.
42 ##
43 ## With no output arguments the results are immediately plotted. With two
44 ## output arguments the 2-D plot data is returned. The data can subsequently
45 ## be plotted manually with @code{plot (@var{x}, @var{y}).
46 ##
47 ## Example:
33 ## 48 ##
34 ## @example 49 ## @example
35 ## @group 50 ## @group
36 ## fplot ("cos", [0, 2*pi]) 51 ## fplot (@cos, [0, 2*pi])
37 ## fplot ("[cos(x), sin(x)]", [0, 2*pi]) 52 ## fplot ("[cos(x), sin(x)]", [0, 2*pi])
38 ## @end group 53 ## @end group
39 ## @end example 54 ## @end example
55 ##
56 ## Note: @code{fplot} works best with continuous functions. Functions with
57 ## discontinuities are unlikely to plot well. This restriction may be
58 ## removed in the future.
40 ## @seealso{plot} 59 ## @seealso{plot}
41 ## @end deftypefn 60 ## @end deftypefn
42 61
43 ## Author: Paul Kienzle <pkienzle@users.sf.net> 62 ## Author: Paul Kienzle <pkienzle@users.sf.net>
44 63
45 function fplot (fn, limits, n = 0.002, fmt = "") 64 function [X, Y] = fplot (fn, limits, varargin)
46 65
47 if (nargin < 2 || nargin > 4) 66 if (nargin < 2 || nargin > 5)
48 print_usage (); 67 print_usage ();
49 endif
50
51 if (iscomplex (limits) || (numel (limits) != 2 && numel (limits) != 4))
52 error ("fplot: LIMITS must be a real vector with 2 or 4 elements");
53 endif
54
55 if (ischar (n))
56 fmt = n;
57 n = 0.002;
58 endif 68 endif
59 69
60 if (strcmp (typeinfo (fn), "inline function")) 70 if (strcmp (typeinfo (fn), "inline function"))
61 fn = vectorize (fn); 71 fn = vectorize (fn);
62 nam = formula (fn); 72 nam = formula (fn);
69 nam = formula (fn); 79 nam = formula (fn);
70 else 80 else
71 error ("fplot: FN must be a function handle, inline function, or string"); 81 error ("fplot: FN must be a function handle, inline function, or string");
72 endif 82 endif
73 83
74 if (n != fix (n)) 84 if (iscomplex (limits) || (numel (limits) != 2 && numel (limits) != 4))
75 tol = n; 85 error ("fplot: LIMITS must be a real vector with 2 or 4 elements");
86 endif
87
88 n = 5;
89 tol = 2e-3;
90 fmt = "";
91 for i = 1:numel (varargin)
92 arg = varargin{i};
93 if (ischar (arg))
94 fmt = arg;
95 elseif (isnumeric (arg) && isscalar (arg) && arg > 0)
96 if (arg == fix (arg))
97 n = arg;
98 else
99 tol = arg;
100 endif
101 else
102 error ("fplot: bad input in position %d", i+2);
103 endif
104 endfor
105
106 if (n != 5)
107 ## n was specified
108 x0 = linspace (limits(1), limits(2), n/2 + 1)';
109 y0 = feval (fn, x0);
110 x = linspace (limits(1), limits(2), n)';
111 y = feval (fn, x);
112 else
76 x0 = linspace (limits(1), limits(2), 5)'; 113 x0 = linspace (limits(1), limits(2), 5)';
77 y0 = feval (fn, x0); 114 y0 = feval (fn, x0);
78 err0 = Inf;
79 n = 8; 115 n = 8;
80 x = linspace (limits(1), limits(2), n)'; 116 x = linspace (limits(1), limits(2), n)';
81 y = feval (fn, x); 117 y = feval (fn, x);
82 118 endif
83 if (! size_equal (x0, y0)) 119
84 ## FN is a constant value function 120 if (rows (x0) != rows (y0))
85 y0 = repmat (y0, size (x0)); 121 ## FN is a constant value function
86 y = repmat (y, size (x)); 122 y0 = repmat (y0, size (x0));
123 y = repmat (y, size (x));
124 endif
125
126 err0 = Inf;
127
128 ## FIXME: This algorithm should really use adaptive scaling as the
129 ## the numerical quadrature algorithms do so that extra points are
130 ## used where they are needed and not spread evenly over the entire
131 ## x-range. Try any function with a discontinuity such as
132 ## fplot (@tan, [-2, 2]) or fplot ("1./x", [-3, 2]) to see the
133 ## problems with the current solution.
134
135 while (n < 2^18) # Something is wrong if we need more than 250K points
136 yi = interp1 (x0, y0, x, "linear");
137 ## relative error calculation using average of [yi,y] as reference
138 ## since neither estimate is known a priori to be better than the other.
139 err = 0.5 * max (abs ((yi - y) ./ (yi + y))(:));
140 if (err < tol || abs (err - err0) < tol/2)
141 ## Either relative tolerance has been met OR
142 ## algorithm has stopped making any reasonable progress per iteration.
143 break;
87 endif 144 endif
88 145 x0 = x;
89 while (n < 2 .^ 20) 146 y0 = y;
90 y00 = interp1 (x0, y0, x, "linear"); 147 err0 = err;
91 err = 0.5 * max (abs ((y00 - y) ./ (y00 + y))(:)); 148 n = 2 * (n - 1) + 1;
92 if (err == err0 || err < tol)
93 break;
94 endif
95 x0 = x;
96 y0 = y;
97 err0 = err;
98 n = 2 * (n - 1) + 1;
99 x = linspace (limits(1), limits(2), n)';
100 y = feval (fn, x);
101 endwhile
102 else
103 x = linspace (limits(1), limits(2), n)'; 149 x = linspace (limits(1), limits(2), n)';
104 y = feval (fn, x); 150 y = feval (fn, x);
105 endif 151 endwhile
106 152
107 plot (x, y, fmt); 153 if (nargout == 2)
108 154 X = x;
109 if (length (limits) > 2) 155 Y = y;
156 else
157 plot (x, y, fmt);
110 axis (limits); 158 axis (limits);
111 endif 159 if (isvector (y))
112 160 legend (nam);
113 if (isvector (y)) 161 else
114 legend (nam); 162 for i = 1:columns (y)
115 else 163 nams{i} = sprintf ("%s(:,%i)", nam, i);
116 for i = 1:columns (y) 164 endfor
117 nams{i} = sprintf ("%s(:,%i)", nam, i); 165 legend (nams{:});
118 endfor 166 endif
119 legend (nams{:});
120 endif 167 endif
121 168
122 endfunction 169 endfunction
123 170
124 171
125 %!demo 172 %!demo
126 %! clf; 173 %! clf;
127 %! fplot ('cos', [0, 2*pi]); 174 %! fplot (@cos, [0, 2*pi]);
128 175
129 %!demo 176 %!demo
130 %! clf; 177 %! clf;
131 %! fplot ('[cos(x), sin(x)]', [0, 2*pi]); 178 %! fplot ('[cos(x), sin(x)]', [0, 2*pi]);
132 179
180 %!demo
181 %! clf;
182 %! ## sinc function
183 %! fh = @(x) sin (x*pi) ./ (x*pi);
184 %! fplot (fh, [-5, 5]);
185
186 %!test
187 %! [x, y] = fplot ("[cos(x), sin(x)]", [0, 2*pi]);
188 %! assert (columns (y) == 2);
189 %! assert (rows (x) == rows (y));
190 %! assert (y, [cos(x), sin(x)], -2e-3);
191
133 %% Test input validation 192 %% Test input validation
134 %!error fplot (1) 193 %!error fplot (1)
135 %!error fplot (1,2,3,4,5) 194 %!error fplot (1,2,3,4,5,6)
195 %!error <FN must be a function handle> fplot (1, [0 1])
136 %!error <LIMITS must be a real vector> fplot (@cos, [i, 2*i]) 196 %!error <LIMITS must be a real vector> fplot (@cos, [i, 2*i])
137 %!error <LIMITS must be a real vector with 2 or 4> fplot (@cos, [1]) 197 %!error <LIMITS must be a real vector with 2 or 4> fplot (@cos, [1])
138 %!error <LIMITS must be a real vector with 2 or 4> fplot (@cos, [1 2 3]) 198 %!error <LIMITS must be a real vector with 2 or 4> fplot (@cos, [1 2 3])
139 %!error <FN must be a function handle> fplot (1, [0 1]) 199 %!error <bad input in position 3> fplot (@cos,[-1,1], {1})
140 200