Borrar filtros
Borrar filtros

Finding the first 50 roots

4 visualizaciones (últimos 30 días)
Jason
Jason el 29 de Mzo. de 2011
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

Respuesta aceptada

Matt Fig
Matt Fig el 29 de Mzo. de 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 comentarios
Jason
Jason el 29 de Mzo. de 2011
Oh and how would I put all the z(n) into a matrix at the end of the computation?
Jason
Jason el 30 de Mzo. de 2011
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.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Loops and Conditional Statements en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by