1
% fzeros_complex.m
2
%
3
% USAGE:
4
%   my_zeros = fzeros_complex(f, z, [optional: print_output])
5
%
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;
8
%     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.
10
%     To find zeros of a function R -> C, set im_min = im_max = 0.
11
% print_output = if this parameter exists, the routine will print some info messages and plot the zeroes found.
12
%
13
% JJ 2011-01-27
14
%
15
function my_zeros = fzeros_complex(f, zig, print_output)
16
	if nargin < 2
17
		error('fzeros_complex.m: parameters required. See "help fzeros_complex" for usage.\n');
18
	end
19
20
	% -------------------------------------------
21
	% parameters
22
	% -------------------------------------------
23
24
    mode = '2d';
25
    if size(zig) == [1,4]
26
		minx = zig(1);
27
		maxx = zig(2);
28
		miny = zig(3);
29
		maxy = zig(4);
30
31
		have_igs = 1;
32
		have_range = 1;
33
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
47
48
%		n_points = 300; % for random initial guess generator
49
	else
50
		have_igs = 1;
51
		have_range = 0;
52
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
62
	if nargin < 3
63
		print_output = 0;
64
	else
65
		print_output = 1;
66
	end
67
68
	max_it = 250;  % max iterations for fminsearch per initial guess
69
	my_tol = 1e-8; % tolerance below which everything is considered zero
70
71
	t_all_start = cputime;
72
73
	% ----------------------------------
74
	% Locate the zeroes
75
	% ----------------------------------
76
77
	% Find the zeroes using fminsearch (Nelder-Mead).
78
	%
79
	t_start = cputime;
80
81
	my_zeros = [];
82
	k = 1;
83
84
	% ||f||^2 = f * conj(f)
85
	normf = @(xy) f(xy(1) + 1i*xy(2)).*conj(f(xy(1) + 1i*xy(2)));
86
87
    % ||f||^2, f : R -> C
88
    normf_real = @(x) f(x).*conj(f(x));
89
90
	% TolX must be really small so that this works even for good initial guesses.
91
	% (See failures, if any, reported by the "candidate" DEBUG message below to find out if the value here is not small enough.)
92
	%
93
	if print_output ~= 0
94
		fprintf(1, 'Trying to find zeroes starting from %d initial guesses...\n', n_points);
95
	end
96
	OPTIONS = optimset('Display', 'off', 'TolFun', my_tol, 'TolX', my_tol * 1e-4, 'MaxIter', max_it);
97
	for j = 1:n_points
98
		if have_igs ~= 0
99
			x = real(zig(j));
100
			y = imag(zig(j));
101
		else
102
			% Generate a random initial guess.
103
			%
104
			x = minx + (maxx - minx) * (rand(1));
105
			y = miny + (maxy - miny) * (rand(1));
106
		end
107
108
		% Minimize ||f||^2 = f * conj(f)
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
116
117
		% DEBUG
118
		if have_igs ~= 0  &&  sv > my_tol  &&  print_output ~= 0
119
			fprintf(1, 'Found local minimum (not accepted as a zero) at z = %0.6g + %0.6g i. f(z) = %0.6g > %0.6g.\n', real(s), imag(s), sv, my_tol);
120
		end
121
122
		% Check if we found a zero.
123
		%
124
		% (The EXITFLAG condition checks that the algorithm found a local minimum,
125
		%  and the second check verifies that it is a zero.)
126
		%
127
		if EXITFLAG == 1  &&  sv <= my_tol
128
			% Aesthetic fix: filter smaller-than-epsilon components out from s.
129
			%
130
			if(abs(real(s)) < my_tol)
131
				s = imag(s) * 1i;
132
			end
133
			if(abs(imag(s)) < my_tol)
134
				s = real(s);
135
			end
136
137
			% Find which so-far-known zero is closest to the one we found
138
			%
139
			[dummy, match_index] = min(abs(my_zeros - s));
140
			
141
			% Always save the first zero.
142
			% Otherwise, save if this zero is farther than my_tol from any known ones.
143
			%
144
			if k == 1  ||  abs(my_zeros(match_index) - s) > my_tol
145
				if print_output ~= 0
146
					fprintf(1, 'Found zero at z = %0.6g + %0.6g i\n', real(s), imag(s));
147
				end
148
				my_zeros(k) = s;
149
				k = k + 1;
150
			end
151
		else
152
%			fprintf(1, 'No zero found for this starting guess (minimum SV = %0.6g)\n', sv); % DEBUG (but shouldn't happen)
153
		end
154
	end
155
156
	if print_output ~= 0
157
		fprintf(1, 'Time taken: %0.3f seconds\n', cputime-t_start);
158
		fprintf(1, 'Found %d zeroes.\n', prod(size(my_zeros)));
159
	end
160
161
	% ----------------------------------
162
	% Phase 3: plot just the zeroes
163
	% ----------------------------------
164
165
	if print_output ~= 0
166
		fprintf(1, 'Plotting zeroes...\n');
167
168
		% plot ||f||^2
169
		if have_range
170
			figure(2);
171
			clf;
172
			nx = 201;
173
			ny = nx;
174
175
			minx = min(minx, min(real(my_zeros)));
176
			maxx = max(maxx, max(real(my_zeros)));
177
			miny = min(miny, min(imag(my_zeros)));
178
			maxy = max(maxy, max(imag(my_zeros)));
179
180
            if strcmp(mode,'1d')
181
    			xx = linspace(minx, maxx, nx);
182
                plot(xx, normf_real(xx), 'k-');
183
184
                title('||f||^2');
185
                xlabel('x');
186
                ylabel('||f||^2');
187
188
                hold on;
189
                plot(my_zeros, zeros(size(my_zeros)), 'ko');
190
                hold off;
191
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)');
205
206
                xx = linspace(minx, maxx, nx);
207
                yy = linspace(miny, maxy, ny);
208
                [X,Y] = meshgrid(xx,yy);
209
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
230
		end
231
	
232
		fprintf(1, 'Done.\n');
233
		fprintf(1, 'Time taken in total: %0.3f seconds\n', cputime-t_all_start);
234
	end
235
end