Finding the first 50 roots
4 views (last 30 days)
Show older comments
The below code was posted by Matt Fig as a response in a thread last fall. I have a similar function that I need to find the first 50 roots. How would you modify the code to accommodate for the change?
function [z] = find_roots()
x = 0.01:0.01:(20*pi/2)-0.01;
plot(x,tan(x)+5.6*x),
grid on
f = @(x) tan(x)+5.6*x;
N = 250; % Fineness of the search, bigger is finer
z = zeros(1,N);
for jj = 1:30
for n=1:N
try
z(n) = fzero(f,jj+[(n-1)/N n/N]);
if abs(f(z(n)))>1e-8
z(n) = 0;
else
fprintf('Found a root at: %.5f\n',z(n))
end
catch
end
end
end
x = 0:pi/50:10*pi;
y = f(x);
plot(z,zeros(1,length(z)),'o',x,y,'-')
line([0 10*pi],[0 0],'color','black')
axis([0 10*pi -1000.0 3000.0])
z = sort(z(z~=0)); % Return only the needed values
0 Comments
Accepted Answer
Matt Fig
on 29 Mar 2011
I think the idea there was to incrementally search for the roots. If you need a certain number of roots, you could use a WHILE loop that evaluates on the root counter. It often helps to know some details of the function in question. Can you post it?
Here is how I would modify the function.
function [z] = find_roots()
f = @(x) x.*cot(x)+99;
NR = 50; % The number of roots to find.
N = 150; % Fineness of the search, bigger is finer z = zeros(1,N);
rcnt = 1; % The root counter
cnt = 0; % The loop counter
while rcnt<=NR
cnt = cnt + 1;
for n = 1:N
try z(rcnt) = fzero(f,cnt+[(n-1)/N n/N]);
if abs(f(z(rcnt)))>1e-8
z(rcnt) = nan; % FZERO was messing with us
else
fprintf('Found a root at: %.5f\n',z(rcnt))
rcnt = rcnt + 1;
break % Put here only if root spacing is >1!
end
catch
end
end
end
x = 0:pi/50:max(z);
plot(x,f(x),'-',z,zeros(1,length(z)),'o')
line([0 max(z)],[0 0],'color','black')
axis([0 max(z) -1000.0 3000.0])
z = z(~isnan(z)); % Return only the needed values
3 Comments
More Answers (0)
See Also
Categories
Find more on Filter Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!