Getting multiple roots of a function with infinite roots

Is there an "easy" way to get multiple roots of a function that has infinite roots? For example solve(sin(x)=0) gives 0, but I also want pi, 2pi, etc.

Do I have to have initial guesses for each solution?

Specifically, I need "an infinite number" (a hundred or so) of positive roots for z*tan(z) = c.

How would I go about doing that?

Answer 1

[Star Strider]

‘Do I have to have initial guesses for each solution?’

For numeric solvers such as fzero, yes. They only look for the closest zero to where you tell them to start looking.

‘How would I go about doing that?’

I’ve used a loop, storing each solution in a vector, then the uniquetol (link) function (in R2015a and later) to eliminate nearly duplicate entries.

[Sentient6]

Problem is, I don't know the roots, but I think they should be more or less close to each other, so the only solution that comes to mind would be really long vector of initial guesses with small steps that will most likely result in a lot of duplicate answers, then using your suggested function. I was hoping there would be something more efficient...

[Star Strider]

To find the approximate real roots, evaluate the function over your region-of-interest, then use this little utility function to get the approximate indices of the roots:

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Approximate Zero-Crossing Indices Of Argument Vector

Then find the x-values associated with those indices, and use those for your starting points.

So if your function is:

y = f(x);
zx = zci(y);
xvals = x(zx);

and go from there.

That’s the best I can do to improve the efficiency of an otherwise random search.

Also, fzero does not do complex roots, while fsolve does, even though you may need to provide it with a complex initial guess.

Answer 2

[David Goodmanson]

Hi Sentient6,

Here is a way that does them all in one go using Newton's method. It's good for either sign of C. x*tanx is an even function so if xn is a root then -xn is a root.

C = 10;
n = 1:50;                       % number of roots
x = n*pi + atan(C./(n*pi));     % initial estimate
for k = 1:20
  xnew = x - (x.*tan(x)-C)./(tan(x)+x./cos(x).^2);
  x = xnew;
end
xn = xnew  

max(abs(xn.*tan(xn)  - C))     % check

That gets them all except for the smallest one when C is positive, which surprisingly seems to be the hardest one to estimate.

C = 10
if C <= 2, x = sqrt(C); else x = atan(C); end
for k = 1:20
  xnew = x - (x.*tan(x)-C)./(tan(x)+x./cos(x).^2);
  x = xnew;
end
x0 = xnew

Answer 3

[shruti hande]

Find the real roots of by plotting the function over the interval