Positive roots of trigonometric equation.
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Risposta accettata
Matt J
il 10 Ott 2014
Modificato: Matt J
il 10 Ott 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 Commenti
John D'Errico
il 10 Ott 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.
Più risposte (1)
Star Strider
il 10 Ott 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 Commenti
Vedere anche
Categorie
Scopri di più su Historical Contests in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!