How to use fzero properly

I've to search for the roots of an equation (the energy levels of multiple quantum well using the propagator method). If I search in all the energy range with a step of 0.01, the run last very long and it doesn't find all the roots. For example in order to find the first 2 energy levels (and having a code that doesn't last too long) I've used this code

if S>=0.75
 for i=0.71544:(0.000001/(N)):0.71546
      try
         liv=fzero(y,[i,i+(0.000001/(N))]);
         LIV(c)=liv;
         c=c+1;
      end
  end 
elseif S>0.5 && S<0.75

    for i=0.7153:(0.00001/(N)):0.7156
        try
           liv=fzero(y,[i,i+(0.00001/(N))]);
           LIV(c)=liv;
           c=c+1;
        end
    end
elseif S>=0.4 && S<=0.5

    for i=0.714:(0.0001/(N)):0.716
        try
           liv=fzero(y,[i,i+(0.0001/(N))]);
           LIV(c)=liv;
           c=c+1;
        end
    end
else
    for i=0.5:(0.01/(N)):0.9
        try
           liv=fzero(y,[i,i+(0.01/(N))]);
           LIV(c)=liv;
           c=c+1;
        end
    end
end

Where S is the distance between the well that is varying and N the number of well. It works but I know where to search for. In another exercise I don't know where to search. How can I scan in all the energy range [0 V0] and using

    for i=0:(0.01/N):(V0)
   try
           liv=fzero(y,[i,i+(0.01/N)]);
          LIV(c)=liv;
          c=c+1;
    end
end

I don't find all the roots

[Image Analyst]

What's the point of the try? Can you post a picture of the function plotted, so we know what we're up against?

[Simo]

the first is the plot using linspace(0,1). and there are no zeros

the second is the same plot using l=0:0.0001:1;

there should be 2 zeros very close (the first two energy levels of a double potential well)

[darova]

Can you attach the original equations?

[Simo]

y=@(E) eq(E,V1);

function [y]=eq(E,V1)
global h_bar_SQ V0 W S N
beta=(sqrt(2*(V0-E)/h_bar_SQ));
alpha=(sqrt(2*E/h_bar_SQ));
Pallowed=[cos(alpha*(W)),(1/alpha)*sin(alpha*(W));-alpha*sin(alpha*(W)),cos(alpha*(W))];
Pforbidden=[cosh(beta*S),(1/beta)*sinh(beta*S);beta*sinh(beta*S),cosh(beta*S)];

alphai=(sqrt(2*(E-(V0-V1))/h_bar_SQ));
Pallowedi=[cos(alphai*(W)),(1/alphai)*sin(alphai*(W));-alphai*sin(alphai*(W)),cos(alphai*(W))];

PSI=[1;beta];


for i=1:2
PSI=Pallowed*PSI; %well
PSI=Pforbidden*PSI;
end

PSI=Pallowedi*PSI; %imp
PSI=Pforbidden*PSI;

for i =1:1
PSI=Pallowed*PSI; %well
PSI=Pforbidden*PSI;
end
PSI=Pallowed*PSI; %well

y=beta*PSI(1)+PSI(2);
end

[John D'Errico]

Sorry it took so long to write the minor treatise on rootfinding in my answer, but that should now be reasonably complete. To go beyond that would really require a course on numerical methods.

Answer 1

[John D'Errico]

Your question is how to use fzero properly. I won't do your assignment, since that is yours to do.

(Oh, see that I've always been careful and created vectorized functions, so they can be evaluated for multiple inputs at once. This will be highly valuable as you learn to use MATLAB.)

fzero wants to have a bracket where a root exists. A bracket means two end points where the function has opposite signs, although an exact zero at one end point would suffice. That rare case would be really unlikely though. Think opposite signs.

Suppose, for example, I wanted to solve the problem of finding some set of roots of a function, perhaps this might be interesting:

fun1 = @(x) x/20 + sin(1./(x + 0.025));
fplot(fun1,[0,20])

where x is a positive number?

The roots surely cannot be found analytically. I've chosen a function that will have a finite number of roots for positive x. But there will be many of them near zero, and they will be very closely packed.

A slightly different case would be

fun2 = @(x) x + tan(x) - 0.5;
fplot(fun2,[0,50])
axis([0 50 -50 50])
yline(0,'r');

Here there are infinitely many roots, and they will not be findable analytically. A nice feature is the roots will be quite regularly spaced, especially for large x. (I recall a question on Answers where I even derived some sort of asymptotic expansion for the set of solutions to a similar problem.) But between each pair of roots, there lies a singularity, where while the function "crosses" zero and it appears a bracket exists, there is never a point where the function actually takes on the value zero. This is because the function is discontinuous at a point, going from +inf to -inf across the discontinuity.

And I can make fun1 even worse that it is. Thus, find the "first" n roots of the function

fun3 = @(x) sin(1./x);
fplot(fun3,[0,1])
yline(0,'r');

for positive x? The problem with fun3 is there are infinitely many such roots less than any positive value you would choose for x. Those roots get more and more closely packed as you approach 0.

There is yet one more rather nasty case to consider for a function, where a root may be nearly impossible to locate. Consider this problem:

delta = 1.e-14;
a = 1000*rand();
fun4 = @(x) 1 - 2*exp(-(x - a).^2/delta);

fun4 has the characteristic that for x anywhere away from a (which I have not told you) the function is almost identically equal to 1.

fun4(17)
ans =
     1

Using double precision arithmetic, essentially ANY value of x will give exactly 1. fzero will NEVER find a zero crossing, even though there are two suc crossings.

a
a =
        673.2

In this case, they will be vey near 673.2, where the function very, very briefly drops down below zero.

format long g
a
a =
          673.195948213502
          
fun4(673.2)
ans =
     1
     
fun4(a)
ans =
    -1
    
fplot(fun4,[673,674])

0

So even told to plot over the interval [673,674], fplot finds nothing at all remarkable, on what appears to be an incredibly boring function.

In order to find that pair of roots, we would effectively need an incredibly fine search interval, at least if we did not have the keys to the city by knowing the value of a in advance.

So using fzero properly in any general case becomes horribly diffiult, because the correct way to use fzero would be different on a case by case basis. I've even considered if I should write a FEX tool for such a problem, but since it would be so easy to cause it to fail, I'm a bit leery.

The classic advice is to use plot FIRST. That is, plot the function, looking for locations where a root would exist, and then throw fzero at each identified root. This works great, except it requires the person to plot, then apply careful thought to what they see, and then call fzero multiple times. While that can work great, it becomes difficult to automate.

Another idea is to simple evaluate the function at a long sequence of points. Then have the computer search for zero crossings. For example, on fun1, we might do this:

fun1 = @(x) x/20 + sin(1./(x + 0.025));
x0 = linspace(0,20,1000);
y0 = fun1(x0);

ind = find(y0(1:end-1).*y0(2:end) <= 0);
brackets = [x0(ind);x0(ind+1)]';

I chose 20 as the upper limit, because I know that no real solution can ever exist for x > 20. (THINK ABOUT IT!)

ind =
     1     2     3     5     7    15

brackets =
                         0          0.02002002002002
          0.02002002002002          0.04004004004004
          0.04004004004004        0.0600600600600601
        0.0800800800800801           0.1001001001001
          0.12012012012012          0.14014014014014
          0.28028028028028           0.3003003003003

I found 6 brackets that "should" contain a solution, because there was a zero crossing identified in that interval. That is, they will contain at least one solution, assuming there was no singularity in the interval. More about that later.

fplot(fun1,brackets(:,1)')

Sadly, even the first bracket has apparently 5 zeros in it. That tells me I used too coarse of an interval when creating x0. The interval spacing was 20/999=0.02002... And of course the problem becomes completely intractable when we might consider what to do for fun3, since ANY finite spacing near zero for x will be too large for that function.

I might also have used a large random sampling of points in the interval for x0. But this tool will fail, unless my sample is made to be impossibly large.

However, the trick usng find as I show it above is a decent one. It will fail on fun2 of course.

fun2 = @(x) x + tan(x) - 0.5;
x0 = linspace(0,100,1000);
y0 = fun2(x0);

ind = find(y0(1:end-1).*y0(2:end) <= 0);
brackets = [x0(ind);x0(ind+1)]';

The problem here is brackets will now contain intervals where an in to inf crossing occurs. Once such event will be in the vicinity of x==pi/2.

brackets(2,:)
ans =
           1.5015015015015           1.6016016016016
           
fplot(fun2,brackets(2,:))

A zero crossing was observed, but no solution will be found there.

In fact, as long as we sample sufficiently finely on fun2, we will find a true zero crossing between each singularity crossing.

>> [xval, fval,exitflag] = fzero(fun2,brackets(1,:))
xval =
         0.247412484885142
fval =
     0
exitflag =
     1
     
>> [xval, fval,exitflag] = fzero(fun2,brackets(2,:))
xval =
           1.5707963267949
fval =
      1.97893796609522e+15
exitflag =
    -5
    
>> [xval, fval,exitflag] = fzero(fun2,brackets(3,:))
xval =
          2.12300090681457
fval =
      2.22044604925031e-16
exitflag =
     1
     
>> [xval, fval,exitflag] = fzero(fun2,brackets(4,:))
xval =
          4.71238898038469
fval =
           510190062007401
exitflag =
    -5

So fzero will return a solution for each case, though beware the non-solutions it did find. A good trick is to test the exit flags, as usually when fzero converges to a non-solution, it is designed to identify them with a different exitflag, signifying that even though it returns something, this might not be a solution you would be happy to see.

How should we best fund the two roots of fun4? Sadly, that can be almost impossible. fzero does not understand that fun4 was created with the explicit purpose of writing a function that would confound fzero.

a
a =
          673.195948213502
>> fun4(a + [-0.01 0 0.01])
ans =
     1    -1     1

If we even deviate from a by the tiniest amount, fun4 always returns exactly 1, to within the abilities of floating point arithmetic. And if that was not enough, I could have made delta significantly smaller yet, turning this subtle deviation that contains two nearly coincident roots into a virtual Dirac delta function.

Much of what I have said here recommends using a bracket to run fzero. In fact, fzero can often survive even when no bracket was provided.

>> fun3 = @(x) sin(1./x);

>> [xval,fval,exitflag] = fzero(fun3,12)
xval =
        -0.318309886183791
fval =
     -1.22464679914735e-16
exitflag =
     1
     
>> [xval,fval,exitflag] = fzero(fun3,1)
xval =
         0.318309886183791
fval =
      1.22464679914735e-16
exitflag =
     1
     
>> [xval,fval,exitflag] = fzero(fun3,.5)
xval =
         0.318309886183791
fval =
     -3.21624529935327e-16
exitflag =
     1
     
>> [xval,fval,exitflag] = fzero(fun3,.25)
xval =
         0.318309886183791
fval =
      1.22464679914735e-16
exitflag =
     1
     
>> [xval,fval,exitflag] = fzero(fun3,.15)
xval =
         0.159154943091895
fval =
     -2.44929359829471e-16
exitflag =
     1

However, as you see, fzero is not able to be picky and know that a diffferent solution is desired on a subsequent call. Lacking a bracket, fzero hunts around until it manages to find a bracket. That search may not be as intelligent as you might like.

Really, the way to correctly use fzero is to first plot your function. My first rule is: ALWAYS PLOT EVERYTHING. Look at what you see. Think about what you see. Then use what you have learned to accomplish your task. Of course, that can be problematic when you want to perform an automated task.

In terms of using fzero, the best general thing you can do is to understand how a root solver works. Understanding how fzero works is even more valuable, of course but that could require an entire class on rootfinding algorithms, since fzero is a code of moderate sophistication. I've already spent far more time on this answer than almost anybody would want to read. But if you want to at least be sem-intelligent in your use of fzero, I would recommend the sort of search I performed in these examples, thus:

1. Sample your function at a sufficiently fine interval to identify any zero crossings.
2. Use fzero on each identified bracket, returning the exitflag.
3. Test the exitflag. If exitflag was not 1, then consider if a root was truly identified.

The possible values for exitflag from fzero are:

1 fzero found a zero X.

-1 Algorithm terminated by output function.

-3 NaN or Inf function value encountered during search for an interval

containing a sign change.

-4 Complex function value encountered during search for an interval

containing a sign change.

-5 fzero may have converged to a singular point.

-6 fzero can not detect a change in sign of the function.

As you should imagine, any exitflag other than 1 is probably a bad thing.

But the above scheme is what I would recommends as the best practice for a good use of fzero. Note that it still requires some intelligence to be applied, as it forces the user to choose a viable sampling interval. You want to sample finely enough that your initial search will not miss any roots. But too fine of an interval just wastes a lot of effort.

For example, on fun2, where the zero crossing interval is roughly pi/2, then anything significantly less than pi/2 for a sampling interval should be sufficient. However, what sampling interval would be neeeded for fun1? The roots get quite close to each other near x==0. And on fun3, we need an infinitely small sampling interval in theory, since ANY finite interval that contains zero also contains infinitely many roots. Likewise, for fun4, I can arbitrarily specify a value for delta small enough to make any sampling scheme fail with probability near 1.

Sadly, I don't think I've given a perfectly useful answer here, because this answer explains more about how to make any root finder fail more than it explains how to make one succeed.