Commit 5ded279e2722608b5b1876581cd88d71d1795bc0

  • Tree SHA1: e7c8a72
  • Parent SHA1: c2160ed (Simple minimization-based zero finder for functions C -> R (requires sudoku_lhs for generating initial guesses))
  • raw diff | raw patch
Support both C->C (2d mode) and R->C (1d mode) functions.
  • fzeros_complex.m 140 --------------------------------------------------++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  • Diff rendering mode:
  • inline
  • side by side

fzeros_complex.m

3% USAGE:3% USAGE:
4% my_zeros = fzeros_complex(f, z, [optional: print_output])4% my_zeros = fzeros_complex(f, z, [optional: print_output])
5%5%
6% f = function handle to a mathematical function C -> R
6% f = function handle to a mathematical function C -> C (or R -> C)
7% z = if a column vector of complex numbers, use them (and nothing else) as initial guesses;7% z = if a column vector of complex numbers, use them (and nothing else) as initial guesses;
8% if a row vector of length 4, use them as [re_min re_max im_min im_max] of the range8% if a row vector of length 4, use them as [re_min re_max im_min im_max] of the range
9% in which initial guesses will be generated.9% in which initial guesses will be generated.
10% To find zeros of a function R -> C, set im_min = im_max = 0.
10% print_output = if this parameter exists, the routine will print some info messages and plot the zeroes found.11% print_output = if this parameter exists, the routine will print some info messages and plot the zeroes found.
11%12%
12% JJ 2011-01-2713% JJ 2011-01-27
21 % parameters21 % parameters
22 % -------------------------------------------22 % -------------------------------------------
2323
24 if size(zig) == [1,4]
24 mode = '2d';
25 if size(zig) == [1,4]
25 minx = zig(1);26 minx = zig(1);
26 maxx = zig(2);27 maxx = zig(2);
27 miny = zig(3);28 miny = zig(3);
31 have_igs = 1;31 have_igs = 1;
32 have_range = 1;32 have_range = 1;
3333
34 k = 20; % boxes for sudoku (per coordinate axis)
35 [zig,m] = sudoku_lhs(2,k,1);
36 n_points = k*m;
37 zig(:,1) = minx + (maxx - minx) * ( (zig(:,1) - 1) / (n_points - 1) );
38 zig(:,2) = miny + (maxy - miny) * ( (zig(:,2) - 1) / (n_points - 1) );
34 % R -> C
35 if miny == 0 && maxy == 0
36 mode = '1d';
37 k = 30;
38 zig = linspace(minx,maxx,k);
39 n_points = k;
40 else % C -> C
41 k = 10; % boxes for sudoku (per coordinate axis)
42 [zig,m] = sudoku_lhs(2,k,1); % N,k,m k*m = n*(k^N)
43 n_points = k*m;
44 zig(:,1) = minx + (maxx - minx) * ( (zig(:,1) - 1) / (n_points - 1) );
45 zig(:,2) = miny + (maxy - miny) * ( (zig(:,2) - 1) / (n_points - 1) );
46 end
3947
40% n_points = 300; % for random initial guess generator48% n_points = 300; % for random initial guess generator
41 else49 else
42 have_igs = 1;50 have_igs = 1;
43 have_range = 0;51 have_range = 0;
44 n_points = length(zig);
45 end
4652
53 m = min(imag(zig(:)));
54 M = max(imag(zig(:)));
55 if m == 0 && M == 0
56 mode = '1d';
57 end
58
59 n_points = length(zig);
60 end
61
47 if nargin < 362 if nargin < 3
48 print_output = 0;63 print_output = 0;
49 else64 else
82 k = 1;82 k = 1;
8383
84 % ||f||^2 = f * conj(f)84 % ||f||^2 = f * conj(f)
85 normf = @(xy) f(xy(1) + i*xy(2)).*conj(f(xy(1) + i*xy(2)));
85 normf = @(xy) f(xy(1) + 1i*xy(2)).*conj(f(xy(1) + 1i*xy(2)));
8686
87 % ||f||^2, f : R -> C
88 normf_real = @(x) f(x).*conj(f(x));
89
87 % TolX must be really small so that this works even for good initial guesses.90 % TolX must be really small so that this works even for good initial guesses.
88 % (See failures, if any, reported by the "candidate" DEBUG message below to find out if the value here is not small enough.)91 % (See failures, if any, reported by the "candidate" DEBUG message below to find out if the value here is not small enough.)
89 %92 %
106 end106 end
107107
108 % Minimize ||f||^2 = f * conj(f)108 % Minimize ||f||^2 = f * conj(f)
109 [S, sv, EXITFLAG] = fminsearch(normf, [x y], OPTIONS);
110 s = S(1) + i * S(2);
109 if strcmp(mode,'1d')
110 [S, sv, EXITFLAG] = fminsearch(normf_real, x, OPTIONS);
111 s = S(1);
112 else
113 [S, sv, EXITFLAG] = fminsearch(normf, [x y], OPTIONS);
114 s = S(1) + 1i * S(2);
115 end
111116
112 % DEBUG117 % DEBUG
113 if have_igs ~= 0 && sv > my_tol && print_output ~= 0118 if have_igs ~= 0 && sv > my_tol && print_output ~= 0
124 % (The EXITFLAG condition checks that the algorithm found a local minimum,124 % (The EXITFLAG condition checks that the algorithm found a local minimum,
125 % and the second check verifies that it is a zero.)125 % and the second check verifies that it is a zero.)
126 %126 %
127 if EXITFLAG == 1 & sv <= my_tol
127 if EXITFLAG == 1 && sv <= my_tol
128 % Aesthetic fix: filter smaller-than-epsilon components out from s.128 % Aesthetic fix: filter smaller-than-epsilon components out from s.
129 %129 %
130 if(abs(real(s)) < my_tol)130 if(abs(real(s)) < my_tol)
131 s = imag(s) * i;
131 s = imag(s) * 1i;
132 end132 end
133 if(abs(imag(s)) < my_tol)133 if(abs(imag(s)) < my_tol)
134 s = real(s);134 s = real(s);
141 % Always save the first zero.141 % Always save the first zero.
142 % Otherwise, save if this zero is farther than my_tol from any known ones.142 % Otherwise, save if this zero is farther than my_tol from any known ones.
143 %143 %
144 if k == 1 | abs(my_zeros(match_index) - s) > my_tol
144 if k == 1 || abs(my_zeros(match_index) - s) > my_tol
145 if print_output ~= 0145 if print_output ~= 0
146 fprintf(1, 'Found zero at z = %0.6g + %0.6g i\n', real(s), imag(s));146 fprintf(1, 'Found zero at z = %0.6g + %0.6g i\n', real(s), imag(s));
147 end147 end
165 if print_output ~= 0165 if print_output ~= 0
166 fprintf(1, 'Plotting zeroes...\n');166 fprintf(1, 'Plotting zeroes...\n');
167167
168 figure(1);
169 clf;
170 plot(real(my_zeros), imag(my_zeros), 'ko');
171 grid on;
172 if have_range
173 axis([minx maxx miny maxy]);
174 end
175 my_title = sprintf('Zeroes of f');
176 title(my_title);
177 xlabel('Re(z)');
178 ylabel('Im(z)');
179
180 % plot ||f||^2168 % plot ||f||^2
181 if have_range169 if have_range
182 figure(2);170 figure(2);
177 miny = min(miny, min(imag(my_zeros)));177 miny = min(miny, min(imag(my_zeros)));
178 maxy = max(maxy, max(imag(my_zeros)));178 maxy = max(maxy, max(imag(my_zeros)));
179179
180 xx = linspace(minx, maxx, nx);
181 yy = linspace(miny, maxy, ny);
182 [X,Y] = meshgrid(xx,yy);
180 if strcmp(mode,'1d')
181 xx = linspace(minx, maxx, nx);
182 plot(xx, normf_real(xx), 'k-');
183183
184 Z = nan(size(X));
185 maxj = prod(size(X));
186 for j = 1:maxj
187 Z(j) = normf([X(j) Y(j)]);
188 end
184 title('||f||^2');
185 xlabel('x');
186 ylabel('||f||^2');
189187
190 surf(X,Y, Z, 'EdgeColor', 'none');
191 view(0,90);
192 title('||f||^2');
193 xlabel('Re(z)');
194 ylabel('Im(z)');
188 hold on;
189 plot(my_zeros, zeros(size(my_zeros)), 'ko');
190 hold off;
195191
196 hold on;
197 plot3(real(my_zeros), imag(my_zeros), 0.1+zeros(size(my_zeros)), 'wo');
198 hold off;
192 axis tight;
193 else
194 figure(1);
195 clf;
196 plot(real(my_zeros), imag(my_zeros), 'ko');
197 grid on;
198 if have_range
199 axis([minx maxx miny maxy]);
200 end
201 my_title = sprintf('Zeroes of f');
202 title(my_title);
203 xlabel('Re(z)');
204 ylabel('Im(z)');
199205
200 colorbar;
206 xx = linspace(minx, maxx, nx);
207 yy = linspace(miny, maxy, ny);
208 [X,Y] = meshgrid(xx,yy);
201209
202 axis tight;
210 Z = nan(size(X));
211 maxj = prod(size(X));
212 for j = 1:maxj
213 Z(j) = normf([X(j) Y(j)]);
214 end
215
216 surf(X,Y, Z, 'EdgeColor', 'none');
217 view(0,90);
218 title('||f||^2');
219 xlabel('Re(z)');
220 ylabel('Im(z)');
221
222 hold on;
223 plot3(real(my_zeros), imag(my_zeros), 0.1+zeros(size(my_zeros)), 'wo');
224 hold off;
225
226 colorbar;
227
228 axis tight;
229 end
203 end230 end
204 231
205 fprintf(1, 'Done.\n');232 fprintf(1, 'Done.\n');