Difference between revisions of "MATLAB:Fzero/Examples"

From PrattWiki
Jump to navigation Jump to search
 
Line 41: Line 41:
 
or use an initial valid bracket.
 
or use an initial valid bracket.
  
== More complicated - multiple possible answers ===
+
== More complicated - multiple possible answers ==
 
The fzero command can only find one root at a time, so you may need to run it several times to calculate all the answers.  As an example, if there is some function:
 
The fzero command can only find one root at a time, so you may need to run it several times to calculate all the answers.  As an example, if there is some function:
 
<center><math>g(y)=(y+1)(y-2)=y^2-y-2</math></center>
 
<center><math>g(y)=(y+1)(y-2)=y^2-y-2</math></center>
Line 65: Line 65:
 
[yVal2, gVal2] = fzero(@(yD) g(yD), [1 4])
 
[yVal2, gVal2] = fzero(@(yD) g(yD), [1 4])
 
</source>
 
</source>
 +
 +
== More complicated - when fzero fails ==
 +
When given an initial guess, fzero actually has to calculate its own bracket.  The way it does this is by bumping the independent variable left and right and determining when the sign of the function at the location to the left has a different sign from the location to the right.  To see this in action, you need to tell fzero to display all iterations of its behavior.  For the following examples, again use:
 +
<center><math>g(y)=(y+1)(y-2)=y^2-y-2</math></center>
 +
To see what happens with an initial guess of 0, use:
 +
<source lang=matlab>
 +
g = @(y) (y+1).*(y-2)
 +
[yVal, gVal] = fzero(@(yD) g(yD), 0, optimset('Display', 'iter'))
 +
</source>
 +
where the optimset command allows yuo to change one or more of the many, many options for fzero (and other commands, too). The result will have two parts.  The first part is finding the bracket:
 +
<source lang=text>
 +
Search for an interval around 0 containing a sign change:
 +
Func-count    a          f(a)            b          f(b)        Procedure
 +
    1              0            -2            0            -2  initial interval
 +
    3      -0.0282843      -1.97092    0.0282843      -2.02748  search
 +
    5          -0.04      -1.9584          0.04      -2.0384  search
 +
    7      -0.0565685      -1.94023    0.0565685      -2.05337  search
 +
    9          -0.08      -1.9136          0.08      -2.0736  search
 +
  11      -0.113137      -1.87406      0.113137      -2.10034  search
 +
  13          -0.16      -1.8144          0.16      -2.1344  search
 +
  15      -0.226274      -1.72253      0.226274      -2.17507  search
 +
  17          -0.32      -1.5776          0.32      -2.2176  search
 +
  19      -0.452548      -1.34265      0.452548      -2.24775  search
 +
  21          -0.64      -0.9504          0.64      -2.2304  search
 +
  23      -0.905097    -0.275703      0.905097      -2.0859  search
 +
  24          -1.28        0.9184      0.905097      -2.0859  search
 +
</source>
 +
a and b are set equal to the initial guess, 0.  Then, a is bumped a little left and f(a) is compared to f(b).  Since the are both negative, b is bumped a little '''right''' and again f(a) and f(b) are checked.  That means MATLAB has run the function three times - once at the original value, once at the new a, and once at the new b.  Each row in the table above, MATLAB bumps a, checks, then bumps be, checks.  In the last line, a finally gets bumped far enough left that f(a) has a positive sign and f(b) continues to have a negative sign.  MATLAB has found an appropriate bracket, and will now proceed as if the function call had been:
 +
<source lang=matlab>
 +
[yVal, gVal] = fzero(@(yD) g(yD), [-1.28  0.905097], optimset('Display', 'iter'))
 +
</source>
 +
 +
The second set of outputs shows you MATLAB's routine for finding a root.  It is a closed method based on interpolation and bisection.  It can find very precise answers very quickly.
 +
 +
=== Bracket fail 1 - Jumping Roots===
 +
the problem occurs when MATLAB's bracketing routine fails to find a bracket.  As you will notice in the table above, MATLAB first tries small bumps then uses progressively larger ones.  If your function has two roots, then, and you guess to far away from where the roots are, the bumping algorithm might bump past two roots and never fine either.  For example, the start of the bracketing table for:
 +
<source lang=matlab>
 +
g = @(y) (y+1).*(y-2)
 +
[yVal, gVal] = fzero(@(yD) g(yD), 30, optimset('Display', 'iter'))
 +
</source>
 +
is
 +
<source lang=text>
 +
Search for an interval around 30 containing a sign change:
 +
Func-count    a          f(a)            b          f(b)        Procedure
 +
    1              30          868            30          868  initial interval
 +
    3        29.1515      818.657      30.8485      918.783  search
 +
    5            28.8        798.64          31.2        940.24  search
 +
    7        28.3029      770.754      31.6971      971.006  search
 +
    9            27.6        732.16          32.4      1015.36  search
 +
  11        26.6059      679.267      33.3941      1079.77  search
 +
  13            25.2        607.84          34.8      1174.24  search
 +
  15        23.2118      513.575      36.7882      1314.59  search
 +
  17            20.4        393.76          39.6      1526.56  search
 +
  19        16.4235      251.309      43.5765      1853.33  search
 +
  21            10.8        103.84          49.2      2369.44  search
 +
  23          2.8471      3.25888      57.1529        3207.3  search
 +
  25            -8.4        76.96          68.4      4608.16  search
 +
  27        -24.3058      613.078      84.3058      7021.16  search
 +
  29          -46.8      2235.04        106.8      11297.4  search
 +
</source>
 +
and you can see that the 24th function call, where a gets bumped from 10.8 to -8.4, is the problem.  The left side of the potential bracket has jumped right past both roots, never to return.  The tail end of the table shows what the problem is:
 +
<source lang=text>
 +
...
 +
2043  -2.84423e+153  8.08962e+306  2.84423e+153  8.08962e+306  search
 +
2045  -4.02234e+153  1.61792e+307  4.02234e+153  1.61792e+307  search
 +
2047  -5.68845e+153  3.23585e+307  5.68845e+153  3.23585e+307  search
 +
2049  -8.04468e+153  6.4717e+307  8.04468e+153  6.4717e+307  search
 +
2051  -1.13769e+154  1.29434e+308  1.13769e+154  1.29434e+308  search
 +
Exiting fzero: aborting search for an interval containing a sign change
 +
    because NaN or Inf function value encountered during search.
 +
(Function value at -1.60894e+154 is Inf.)
 +
Check function or try again with a different starting value.
 +
</source>
 +
The key then is to either make sure you start close enough to the root or provide a good initial bracket to begin with.

Revision as of 21:57, 1 October 2012

The following will be a series of more examples using the fzerocommand:

Very Simple

Assume there is some function

\(f(x)=x+1\)

and you are trying to determine the value of \(x\) for which \(f(x)=0\). There are many different ways to use MATLAB and fzero to find the solution:

Put the calculation in the fzero command

If the function only requires one line of code to write, you can put the anonymous function that does the calculation directly in the fzero command. fzero will know it has control only over the dummy variable given in the argument list after the @ symbol. You will also need to give either an initial guess, for example\[x=2\]

[xVal, fVal] = fzero(@(xD) xD+1, 2)

or a valid initial bracket (valid meaning the sign of the function changes from one side to the other):

[xVal, fVal] = fzero(@(xD) xD+1, [-2 2])

If you give an invalid initial bracket, such as

[xVal, fVal] = fzero(@(xD) xD+1, [0 5])

where in this case the values of the function are positive at both \(f(0)\) and \(f(5)\), fzero will not run. Note in the remaining examples below that the initial guess case is generally used; in all those cases, you can substitute in a valid initial bracket.

Create a named anonymous function and then call it

You can also create the anonymous function in advance and then call it in the fzero command. This is useful if you need to use the function again later, say, for a plot:

MyFun = @(x) x+1
[xVal, fVal] = fzero(@(xD) MyFun(xD), 0)

or, again, you can use an initial bracket as long as it is valid.There's extra code here, but again, if you need to use the function for more than just finding zeros, this is convenient.

Put the function in a function file and then call it

If you have a calculation that requires more than one line of code, you can create your own function and then have fzero call it. For example, you could create a file called MyFunFile.m that contains

function out = MyFunFile(in)
    out = in + 1;

abd then, in the command window or in your script, you can write <course lang=matlab> [xVal, fVal] = fzero(@xD) MyFunFile(xD), 0) </course> or use an initial valid bracket.

More complicated - multiple possible answers

The fzero command can only find one root at a time, so you may need to run it several times to calculate all the answers. As an example, if there is some function:

\(g(y)=(y+1)(y-2)=y^2-y-2\)

and you are trying to determine the value(s) of \(y\) for which \(g(y)=0\) you will need to run fzero twice. You can try to give initial guesses close to the root you are looking for:

g = @(y) (y+1).*(y-2)
[yVal1, gVal1] = fzero(@(yD) g(yD), 0)
[yVal2, gVal2] = fzero(@(yD) g(yD), 3)

ends up working. Typically, however, you will want to graph the equation, look for the sign changes, and then use a valid bracket to make sure you are getting the right root. That means first writing something on the order of:

g = @(y) (y+1).*(y-2)
y = linspace(-4, 4, 100);
plot(y, g(y), 'k-')

to see the values or

plot(y, sign(g(y)), 'k-')

to see the sign changes followed by fzero commands that bracket the individual roots. The actual independent values for the bracket may vary; here are some examples:

[yVal1, gVal1] = fzero(@(yD) g(yD), [-2 0])
[yVal2, gVal2] = fzero(@(yD) g(yD), [1 4])

More complicated - when fzero fails

When given an initial guess, fzero actually has to calculate its own bracket. The way it does this is by bumping the independent variable left and right and determining when the sign of the function at the location to the left has a different sign from the location to the right. To see this in action, you need to tell fzero to display all iterations of its behavior. For the following examples, again use:

\(g(y)=(y+1)(y-2)=y^2-y-2\)

To see what happens with an initial guess of 0, use:

g = @(y) (y+1).*(y-2)
[yVal, gVal] = fzero(@(yD) g(yD), 0, optimset('Display', 'iter'))

where the optimset command allows yuo to change one or more of the many, many options for fzero (and other commands, too). The result will have two parts. The first part is finding the bracket:

Search for an interval around 0 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1               0            -2             0            -2   initial interval
    3      -0.0282843      -1.97092     0.0282843      -2.02748   search
    5           -0.04       -1.9584          0.04       -2.0384   search
    7      -0.0565685      -1.94023     0.0565685      -2.05337   search
    9           -0.08       -1.9136          0.08       -2.0736   search
   11       -0.113137      -1.87406      0.113137      -2.10034   search
   13           -0.16       -1.8144          0.16       -2.1344   search
   15       -0.226274      -1.72253      0.226274      -2.17507   search
   17           -0.32       -1.5776          0.32       -2.2176   search
   19       -0.452548      -1.34265      0.452548      -2.24775   search
   21           -0.64       -0.9504          0.64       -2.2304   search
   23       -0.905097     -0.275703      0.905097       -2.0859   search
   24           -1.28        0.9184      0.905097       -2.0859   search

a and b are set equal to the initial guess, 0. Then, a is bumped a little left and f(a) is compared to f(b). Since the are both negative, b is bumped a little right and again f(a) and f(b) are checked. That means MATLAB has run the function three times - once at the original value, once at the new a, and once at the new b. Each row in the table above, MATLAB bumps a, checks, then bumps be, checks. In the last line, a finally gets bumped far enough left that f(a) has a positive sign and f(b) continues to have a negative sign. MATLAB has found an appropriate bracket, and will now proceed as if the function call had been:

[yVal, gVal] = fzero(@(yD) g(yD), [-1.28  0.905097], optimset('Display', 'iter'))

The second set of outputs shows you MATLAB's routine for finding a root. It is a closed method based on interpolation and bisection. It can find very precise answers very quickly.

Bracket fail 1 - Jumping Roots

the problem occurs when MATLAB's bracketing routine fails to find a bracket. As you will notice in the table above, MATLAB first tries small bumps then uses progressively larger ones. If your function has two roots, then, and you guess to far away from where the roots are, the bumping algorithm might bump past two roots and never fine either. For example, the start of the bracketing table for:

g = @(y) (y+1).*(y-2)
[yVal, gVal] = fzero(@(yD) g(yD), 30, optimset('Display', 'iter'))

is

Search for an interval around 30 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1              30           868            30           868   initial interval
    3         29.1515       818.657       30.8485       918.783   search
    5            28.8        798.64          31.2        940.24   search
    7         28.3029       770.754       31.6971       971.006   search
    9            27.6        732.16          32.4       1015.36   search
   11         26.6059       679.267       33.3941       1079.77   search
   13            25.2        607.84          34.8       1174.24   search
   15         23.2118       513.575       36.7882       1314.59   search
   17            20.4        393.76          39.6       1526.56   search
   19         16.4235       251.309       43.5765       1853.33   search
   21            10.8        103.84          49.2       2369.44   search
   23          2.8471       3.25888       57.1529        3207.3   search
   25            -8.4         76.96          68.4       4608.16   search
   27        -24.3058       613.078       84.3058       7021.16   search
   29           -46.8       2235.04         106.8       11297.4   search

and you can see that the 24th function call, where a gets bumped from 10.8 to -8.4, is the problem. The left side of the potential bracket has jumped right past both roots, never to return. The tail end of the table shows what the problem is:

...
 2043   -2.84423e+153  8.08962e+306  2.84423e+153  8.08962e+306   search
 2045   -4.02234e+153  1.61792e+307  4.02234e+153  1.61792e+307   search
 2047   -5.68845e+153  3.23585e+307  5.68845e+153  3.23585e+307   search
 2049   -8.04468e+153   6.4717e+307  8.04468e+153   6.4717e+307   search
 2051   -1.13769e+154  1.29434e+308  1.13769e+154  1.29434e+308   search
Exiting fzero: aborting search for an interval containing a sign change
    because NaN or Inf function value encountered during search.
(Function value at -1.60894e+154 is Inf.)
Check function or try again with a different starting value.

The key then is to either make sure you start close enough to the root or provide a good initial bracket to begin with.