Finding the first 50 roots

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

 채택된 답변

Matt Fig
Matt Fig 2011년 3월 29일

1 개 추천

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

Jason
Jason 2011년 3월 29일
the function is x.*cot(x)+99 = 0
so would you put a counter outside the for loop as in while counter<50
Jason
Jason 2011년 3월 29일
Oh and how would I put all the z(n) into a matrix at the end of the computation?
Jason
Jason 2011년 3월 30일
That is perfect Matt, you are the magic man. Thank you very much. You just saved me several hours of banging my head against a wall.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

도움말 센터File Exchange에서 Marine and Underwater Vehicles에 대해 자세히 알아보기

태그

질문:

2011년 3월 29일

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by