Find first n roots of transcendental equation

I'm trying to find the first n roots of the equation

,

where α is a constant I have already calculated.

I know that the roots are close are essentially equal to the odd numbers, but I need to show it. I have the following code:

a = 0.9092;

myfun = @(lam,a) tan(lam*pi/2)-lam*pi/2+(4/a^2)*(lam*pi/2)^3;   % parameterized function
fun = @(lam) myfun(lam,a);                                      % function of lam alone

lam_i = zeros(100,1);
for n = 1:2:length(lam_i)
    lam_i(n) = fzero(fun,n);
end
lam_i

lam_i = 100×1

1.0000 0 -0.0000 0 0.0000 0 -0.0000 0 -0.0000 0

The problem is that it misses out some of the roots I know it has, i.e. it at one point finds 71, 63, 79, 83 when I want 63,65,67 etc. It also includes the trivial solution numerous times.

Is there a better way to solve an equation like this?

Thanks :)

Answer 1

[John D'Errico]

First, the function you have written does not seem to have multiple roots. In fact, the only solution it seems to have is zero. Here, I have copied and pasted your code, then plotted the resulting function:

a = 0.9092;
myfun = @(lam,a) tan(lam*pi/2)-lam*pi/2+(4/a^2)*(lam*pi/2).^3;   % parameterized function
fun = @(lam) myfun(lam,a);                                      % function of lam alone
fplot(fun,[-20,20])

You can see it crosses the axis only at lam == 0.

But if you had a problem with multiple roots...

Is there a way that will insure you ALWAYS find all roots of a fully general problem with multiple roots on some interval? NO.

In fact, if the problem is not limited to be sufficiently well behaved, thus has limits on the slopes and derivatives of the function, then it can be proved that no numerical method can solve such a problem and not miss any root.

So what does this tell you? That the answer to your question is, sorry, but no.

What you did was to run a numerical method, giving it starting values at integer increments, hoping that would suffice. But it won't, at least not for the case in point. Can you do better? Perhaps.

First, do not give fzero a SINGLE starting value. But give it a pair of values where the function has distinct signs. Then fzero will alwsys find a root between that pair of function values. Is this sufficient to find all roots? No. Because if there are three roots between any bracket that changes sign over that interval, then fzero will find only 1 of those roots.

So a better code would look at a fairly fine discretization. Evaluate the function at EVERY point in that list. If the function changes sign between any consecutive pair, of point, then you have found a bracket to call fzero on. Will this insure you find EVERY root? No. Re-read what I said above. Nothing you do will always succeed. You can only improve the odds that you will succeed, by making a sufficiently fine search.

In the case where you think the root lies roughly at the odd integers, then a search might look like that below:

I'll use an example case where I know there are multiple roots. Here it will be tan(x) == x.

ntarget = 10; % number of roots you want
rootlist = NaN(1,ntarget);
nfound = 0;
x0 = 0;

fun = @(x) tan(x) - x; 

dx = 1/(3*exp(1)); % some small interval, chosen to be not evenly divisible into 1.

while nfound < ntarget
    if sign(fun(x0)) ~= sign(fun(x0+dx))
        nfound = nfound + 1;
        rootlist(nfound) = fzero(fun,x0 + [0,dx]);
    end
    x0 = x0 + dx;
end

rootlist

rootlist = 1×10

0 1.5708 4.4934 4.7124 7.7253 7.8540 10.9041 10.9956 14.0662 14.1372

In the case of tan(x)==x, the roots will tend to lie at increments that differ by approximately pi, as x grows larger. And you can see it does indeed do that. However, you will need to be careful, since fzero also has found the zero crossings at every point pi/2+n*pi. That is not a root, but a singularity of the tangent function, where tan(x) goes from +inf to -inf.

rootlist/pi

ans = 1×10

0 0.5000 1.4303 1.5000 2.4590 2.5000 3.4709 3.5000 4.4774 4.5000

This is a problem that your code will also need to deal with, since any function that has tan in it like that will also have singularities.

Answer 2

[Star Strider]

First, plot the funciton —

a = 0.9092;
myfun = @(lam,a) tan(lam*pi/2)-lam*pi/2+(4/a^2).*(lam*pi/2).^3; % parameterized function
fun = @(lam) myfun(lam,a);                                      % function of lam alone

lamv = linspace(0, 100, 1E+3);

figure
subplot(2,1,1)
plot(lamv, fun(lamv))
grid
subplot(2,1,2)
plot(lamv, fun(lamv))
grid
ylim([-1 1])

There only appear to be two roots, at about 0.9 and about 1.1. (The function value near 0 is a trivial solution.) The second subplot only restricts the y-limits, so any other zero-crossings would be plotted. There are none.

.