find roots

Hi everyone.

I'm struggling with a pretty easy problem (I'm using Matlab since short!): I have to find the (1000 first) roots of the following function:

y = x*cot(x) + L -1

where c is a constant (in my case, usually between 1 and 100). I had to use a few tricks to avoid the undesired solutions that the fzero-function leads to. In the end, my code is the following (I called x beta):

function [beta, beta_null] = f_find_beta(L)

x0 = 0.001;
L = 10;
beta = zeros(1,1000);
beta_null = zeros(1,999);
%betaq = 0.1;

for q = 1:1:1000

      x = linspace( x0, x0 + pi );
      y = x .* cot(x) + L - 1;
      for p = 1:1:99
          if y( p+1 ) ./ y( p ) < 0
              z = fzero( @(x) x .* cot(x) + L - 1, x(p) );
              if abs( (z ./ pi) - (round(z ./ pi)) ) > 0.01 %&& (z - betaq) > 2.0
                  beta(q) = z;
                  %betaq = z;
              else end
          else end
      end
      x0 = x0 + pi ;
end

for n = 1:1:999
    beta_null(n) = beta(n+1) - beta(n);
end

end

beta_null is just a way for me to check my results more quickly. If you plot this vector as a function of its indices, you should get an almost straight line around pi.

It turns out that only the first 34-35 first roots can be calculated with this method without mistake. Then Matlab mistakes a "jump" of the function from +Inf to -Inf with a possible location of a root. I avoided the first mistakes with one of the if-loops (using a property of the x*cot(x)-function: its asymptotes are periodic (but not its roots) ) but the second trick did not improve my .m-file any further (see the "betaq" lines, with the % at the beginning since this method does not work).

Do you know how I could improve this file? In case this file seems to you to be a disaster, do you have another solution? (I also tried using the "findallzeros"-file that one can find really easily by typing "find roots matlab" or "find zeros matlab" in Google, but it does not work either).

Thanks for your help.

Tibo

Answer 1

[Andrew Newell]

The key to this problem is that there is one root in each interval [n*pi,(n+1)*pi]. I'm going to assume that L > 1, otherwise the region near the origin is a special case.

function beta = find_beta(L,nRoots)

f = @(x) x*cot(x)+L-1;

beta = zeros(nRoots,1);

% A small offset is needed to avoid the asymptotes.
smallOffset = 1e-12;

for i=1:nRoots
    beta(i) = fzero(f,[i*pi+smallOffset (i+1)*pi-smallOffset]);
end

[Thibaut]

Damn, that's beautifully simple.

Thanks a lot.

[Andrew Newell]

Thibaut added:

I noticed that the maximal number of roots I can compute with it is restricted by the parameter SmallOffset. With the code you wrote, only 5220 values can be computed. I am now studying a case in which I need to compute 1E4 values of beta. It works with SmallOffset = 1E-11 (instead of 1E-12). Can you explain me why?

It has nothing to do with the exclusion of a root due to the offset, because the roots - in this rage - lie right in the middle of the interval given in the fzero-function (which is also why the code still works with SmallOffset = 1E-11).

[Andrew Newell]

It's rounding error. The floating-point accuracy (which you can see by typing *eps*) is about 2e-16. Now 1e-12 / 5220 / pi is 6e-17, so it's not surprising that MATLAB starts to have trouble distinguishing between just below 5220*pi and just above.

[Thibaut]

Once again, very helpful. Thank you very much.