Finding the first 50 roots

4 次查看(过去 30 天)
Jason
Jason 2011-3-29
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
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
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 个)

类别

Help CenterFile Exchange 中查找有关 Marine and Underwater Vehicles 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by