Difference between revisions of "MATLAB:Fzero/Examples"
m |
|||
| Line 35: | Line 35: | ||
out = in + 1; | out = in + 1; | ||
</source> | </source> | ||
| − | + | and then, in the command window or in your script, you can write | |
<source lang=matlab> | <source lang=matlab> | ||
[xVal, fVal] = fzero(@(xD) MyFunFile(xD), 0) | [xVal, fVal] = fzero(@(xD) MyFunFile(xD), 0) | ||
Revision as of 15:50, 19 October 2014
The following will be a series of more examples using the fzero command:
Contents
Very Simple
Assume there is some function
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:
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:
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.
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 colve
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.
