Positive roots of trigonometric equation.
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Cedric
el 10 de Oct. de 2014
Comentada: Star Strider
el 11 de Oct. de 2014
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.
0 comentarios
Respuesta aceptada
Matt J
el 10 de Oct. de 2014
Editada: Matt J
el 10 de Oct. de 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 comentarios
John D'Errico
el 10 de Oct. de 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.
Más respuestas (1)
Star Strider
el 10 de Oct. de 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 comentarios
Ver también
Categorías
Más información sobre Historical Contests en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!