Positive roots of trigonometric equation.

5 views (last 30 days)
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.

Accepted Answer

Matt J
Matt J on 10 Oct 2014
Edited: Matt J on 10 Oct 2014
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.
  2 Comments
John D'Errico
John D'Errico on 10 Oct 2014
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
Cedric on 11 Oct 2014
Edited: Cedric on 11 Oct 2014
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

Sign in to comment.

More Answers (1)

Star Strider
Star Strider on 10 Oct 2014
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.
  3 Comments
Cedric
Cedric on 11 Oct 2014
Hi Star, thank you for your answer!

Sign in to comment.

Categories

Find more on Historical Contests in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!