Difference between revisions of "MATLAB:Fzero/Examples"

From PrattWiki
Jump to navigation Jump to search
m
 
(6 intermediate revisions by 2 users not shown)
Line 1: Line 1:
−
The following will be a series of more examples using the <code>fzero</code>command:
+
The following will be a series of more examples using the <code>[[MATLAB:Fzero|fzero]]</code> command:
  
 
== Very Simple ==
 
== Very Simple ==
 
Assume there is some function
 
Assume there is some function
−
<center><math>f(x)=x+1</math></center>
+
<center><math>f(x)=x+1\!</math></center>
 
and you are trying to determine the value of <math>x</math> for which <math>f(x)=0</math>.  There are many different ways to use MATLAB and fzero to find the solution:
 
and you are trying to determine the value of <math>x</math> for which <math>f(x)=0</math>.  There are many different ways to use MATLAB and fzero to find the solution:
 
=== Put the calculation in the fzero command ===
 
=== Put the calculation in the fzero command ===
Line 35: Line 35:
 
     out = in + 1;
 
     out = in + 1;
 
</source>
 
</source>
−
abd then, in the command window or in your script, you can write
+
and then, in the command window or in your script, you can write
−
<course lang=matlab>
+
<source lang=matlab>
−
[xVal, fVal] = fzero(@xD) MyFunFile(xD), 0)
+
[xVal, fVal] = fzero(@(xD) MyFunFile(xD), 0)
−
</course>
+
</source>
 
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>
 
and you are trying to determine the value(s) of <math>y</math> for which <math>g(y)=0</math> you will need to run fzero twice.  You can try to give initial guesses close to the root you are looking for:
 
and you are trying to determine the value(s) of <math>y</math> for which <math>g(y)=0</math> you will need to run fzero twice.  You can try to give initial guesses close to the root you are looking for:
 
<source lang=matlab>
 
<source lang=matlab>
Line 68: Line 68:
 
== More complicated - when fzero fails ==
 
== 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:
 
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>
+
<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:
 
To see what happens with an initial guess of 0, use:
 
<source lang=matlab>
 
<source lang=matlab>
Line 92: Line 92:
 
   24          -1.28        0.9184      0.905097      -2.0859  search
 
   24          -1.28        0.9184      0.905097      -2.0859  search
 
</source>
 
</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:
+
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 they 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 b, 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>
 
<source lang=matlab>
 
[yVal, gVal] = fzero(@(yD) g(yD), [-1.28  0.905097], optimset('Display', 'iter'))
 
[yVal, gVal] = fzero(@(yD) g(yD), [-1.28  0.905097], optimset('Display', 'iter'))
Line 139: Line 139:
 
</source>
 
</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.
 
The key then is to either make sure you start close enough to the root or provide a good initial bracket to begin with.
 +
=== Bracket fail 2 - going complex ===
 +
The other way MATLAB's bracketing can fail is if bumping puts the function into a region where the solution becomes complex or imaginary.  Imagine you want to solve
 +
<center><math>h(t)=\sqrt{t}</math></center>
 +
for values of <math>t</math> where <math>h(t)=2</math>.  To begin, you need to change this to a root-finding problem, so you will actually be solving for <math>h(t)-2=0</math>.  You can do this a variety of ways, including:
 +
<source lang=matlab>
 +
h = @(t) sqrt(t)
 +
[tVal, funVal] = fzero(@(tD) h(tD)-2, [0 10])
 +
</source>
 +
and that will work fine.  But if you try:
 +
<source lang=matlab>
 +
h = @(t) sqrt(t)
 +
[tVal, funVal] = fzero(@(tD) h(tD)-2, 1)
 +
</source>
 +
you get an error.  Displaying iterations might help discover the problem:
 +
<source lang=matlab>
 +
h = @(t) sqrt(t)
 +
[tVal, funVal] = fzero(@(tD) h(tD)-2, 1, optimset('Display', 'iter'))
 +
</source>
 +
In this case, the bracketing routine shows:
 +
<source lang=text>
 +
Search for an interval around 1 containing a sign change:
 +
Func-count    a          f(a)            b          f(b)        Procedure
 +
    1              1            -1            1            -1  initial interval
 +
    3        0.971716      -1.01424      1.02828    -0.985956  search
 +
    5            0.96      -1.0202          1.04    -0.980196  search
 +
    7        0.943431      -1.0287      1.05657    -0.972105  search
 +
    9            0.92      -1.04083          1.08      -0.96077  search
 +
  11        0.886863      -1.05827      1.11314    -0.944947  search
 +
  13            0.84      -1.08348          1.16    -0.922967  search
 +
  15        0.773726      -1.12038      1.22627    -0.892627  search
 +
  17            0.68      -1.17538          1.32    -0.851087  search
 +
  19        0.547452      -1.2601      1.45255    -0.794783  search
 +
  21            0.36          -1.4          1.64    -0.719375  search
 +
  23      0.0949033      -1.69194        1.9051    -0.619748  search
 +
Exiting fzero: aborting search for an interval containing a sign change
 +
    because complex function value encountered during search.
 +
(Function value at -0.28 is -2+0.52915i.)
 +
Check function or try again with a different starting value.
 +
</source>
 +
It turns out that the routine bumped a to -0.28, and <math>h(-0.28)-2</math> is a complex number.  MATLAB's fzero routine cannot handle complex roots and thus ends immediately.  This is a case where a good bracket will protect you - making sure the left side of the bracket does not go past 0 will ensure the values are real.

Latest revision as of 15:55, 19 October 2014

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

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;

and then, in the command window or in your script, you can write 

[xVal, fVal] = fzero(@(xD) MyFunFile(xD), 0)

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 they 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 b, 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.

Bracket fail 2 - going complex

The other way MATLAB's bracketing can fail is if bumping puts the function into a region where the solution becomes complex or imaginary. Imagine you want to solve

\(h(t)=\sqrt{t}\)

for values of \(t\) where \(h(t)=2\). To begin, you need to change this to a root-finding problem, so you will actually be solving for \(h(t)-2=0\). You can do this a variety of ways, including: 

h = @(t) sqrt(t)
[tVal, funVal] = fzero(@(tD) h(tD)-2, [0 10])

and that will work fine. But if you try: 

h = @(t) sqrt(t)
[tVal, funVal] = fzero(@(tD) h(tD)-2, 1)

you get an error. Displaying iterations might help discover the problem: 

h = @(t) sqrt(t)
[tVal, funVal] = fzero(@(tD) h(tD)-2, 1, optimset('Display', 'iter'))

In this case, the bracketing routine shows: 

Search for an interval around 1 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1               1            -1             1            -1   initial interval
    3        0.971716      -1.01424       1.02828     -0.985956   search
    5            0.96       -1.0202          1.04     -0.980196   search
    7        0.943431       -1.0287       1.05657     -0.972105   search
    9            0.92      -1.04083          1.08      -0.96077   search
   11        0.886863      -1.05827       1.11314     -0.944947   search
   13            0.84      -1.08348          1.16     -0.922967   search
   15        0.773726      -1.12038       1.22627     -0.892627   search
   17            0.68      -1.17538          1.32     -0.851087   search
   19        0.547452       -1.2601       1.45255     -0.794783   search
   21            0.36          -1.4          1.64     -0.719375   search
   23       0.0949033      -1.69194        1.9051     -0.619748   search
Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at -0.28 is -2+0.52915i.)
Check function or try again with a different starting value.

It turns out that the routine bumped a to -0.28, and \(h(-0.28)-2\) is a complex number. MATLAB's fzero routine cannot handle complex roots and thus ends immediately. This is a case where a good bracket will protect you - making sure the left side of the bracket does not go past 0 will ensure the values are real.

Retrieved from "https://pundit.pratt.duke.edu/piki/index.php?title=MATLAB:Fzero/Examples&oldid=17870"