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 Center および File ExchangeFilter Analysis についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by