Positive roots of trigonometric equation.

Dear all,

What would be a sound and reliable way for finding the first n positive roots of the following?

 a   = 918 ;
 b   = 47 ; 
 c   = 3e-2 ;
 fun = @(x) x .* tan(x) - (a-x.^2) ./ (b + (a-x.^2)*c) ;

My first idea was to use FZERO from a large number of starting values, keep unique (with some tolerance) positive roots, and eliminate those close (k+1/2)pi .. but there has to be a sounder and more reliable approach (?)

Thank you,

Cedric

PS/EDITs:

• And eliminate as well roots of the denominator.

Answer 1

[Matt J]

The function

     f(x) = x .* tan(x)

is continuous and strictly monotonically increasing in every interval [-pi/2,pi/2]+k*pi while the function

     g(x) = (a-x.^2) ./ (b + (a-x.^2)*c)

is continuous and strictly monotonically decreasing everywhere except at its unique pole xp.

Because f and g have opposite monotonicity in each interval [-pi/2,pi/2]+k*pi not containing xp, and because f(x) ranges from -inf to +inf there, the equation

                    f(x)=g(x)

will have exactly one solution in each of those intervals. These solutions are identically the roots of f(x)-g(x). The anomalous interval [-pi/2,pi/2]+k0*pi containing the pole xp can be handled by splitting it into 2 subintervals to the left and right of xp. Each of those subintervals must have exactly one root by the same monotonicity arguments.

Thus, you need simply loop over the first successive n of these intervals,subdividing the k0-th interval appropriately. In each interval, use fzero to find the unique root there.

[John D'Errico]

While I like this answer very much (and up-voted it), I'd want to be more convincing that g(x) is indeed a decreasing function on both sides of the pole

xp = sqrt(22362)/3 = 49.846430831772365391276755808461

As long as one is only interested in positive roots, that is indeed true.

[Cedric]

Thank you Matt for your answer! As parameters can vary (which I didn't mention; and I don't know yet their ranges), I'll have to check out the monotonicity, but it is a small price to pay for a great approach!

Edit:

 sign(dg/dx) = -sign(b) for x >= 0

Answer 2

[Star Strider]

My approach is strictly numeric:

x = linspace(0,20*pi,1E4);              % Define Domain
y = fun(x);                             % Evaluate Range

ycs = y.*circshift(y, [0 -1]);          % Shift by 1 & Multiply
yzx = find(ycs <= 0);                   % Approx Zero Crossing Indices
yzx = yzx(1:2:end);                     % Choose Alternates

ixr = 10;                               % Interpolation Range
for k1 = 1:length(yzx)
    b = [ones(1+ixr*2,1)  x(yzx(k1)-ixr:yzx(k1)+ixr)']\y(yzx(k1)-ixr:yzx(k1)+ixr)';
    rt(k1) = -b(1)/b(2);
end

figure(1)
plot(x, y, '-b')
hold on
% plot(x(yzx), zeros(size(yzx)), 'xm')
plot(rt, zeros(size(rt)), 'pg')
hold off
grid
axis([0  10    [-1  1]*100 ])

This is a simple linear interpolation method that is likely faster than interp1. The number of roots it finds is determined by the second argument of linspace. When I timed it from before linspace to the end of the loop, it took 0.01 seconds on my less-than-frisky laptop to find 21 roots.

[Matt J]

The number of roots it finds is determined by the second argument of linspace.

And the 3rd argument, of course. Had you done

   linspace(0,20*pi,2)

you could find at most 1 root. You also aren't guaranteeing that you will find sequential roots, as far as I can tell, i.e., some roots can be skipped over. You probably could prove it doesn't happen for this particular problem, but that would require more than strictly numeric work.

[Cedric]

Hi Star, thank you for your answer!

[Star Strider]

My pleasure!