A Recipe for Finding Roots of
Finding the value of that satisfies
requires an iterative guess-and-correct approach. Generally this is called a root-finding problem. The standard solution begins rewriting the equation as
and then using a numerical method to find the value of that makes .
Note that we do not use a "solve" method. Novices confuse the solve
function in the symbolic math toolbox with the numerical procedure
of finding a value of that comes close, usually very close,
to giving .
For numerical root-finding, do not use solve
.
In this example, I walk you through the following steps
- Plot to see where it crosses the axis.
- Identify a good bracket interval that contains the root.
- Use bisection, or some other root-finding methods, to refine the guess at the root within the interval.
- Verify that a root was found.
1. Plot
Before we can plot it is helpful to write a simple m-file function to evaluate for any reasonable input, . This m-file function will be needed later in the root-finding process. Here is the code
function z = myfx(x)
z = x - x.^(1/3) - 2;
end
This trivially simple m-file has several important features that are necessary to streamline the root-finding process later
- The code is vectorized, meaning that the input can be a scalar value or a vector.
- The code does only one thing: evaluate given
- The code has a single input
x
and it returns a single output.
With the myfx
function defined in the myfx.m
file, the following statements
create a plot of .
>> x = linspace(0,5);
>> f = myfx(x);
>> plot(x,f,'-')
>> grid('on')
Running the preceding code snippet produces this plot.
Looking at the plot we can see that crosses the axis somewhere between and . In other words, the solution of is in the interval .
2. Identify a Bracket that Contains a Root
The plot of gives us the information we need to initiate an automatic
root-finding procedure, for example, using the bisect
function that implements
the bisection algorithm. Though it is strictly not necessary, we can introduce
some automation into this process with the brackPlot
helper function.
Running a single statement
>> xb = brackplot('myfx',0,5)
xb =
3.4211 3.6842
produces the following plot and stores the bracket interval in the vector, xb
.
Note: Using brackPlot
means that we do not generate a vector of
values with x = linspace(...)
and then manually plotting with
f =myfx(x); plot(x,f)
. brackPlot
does that for us.
3. Refine the Guess at the Root
With a bracket to the root stored in xb
, another one-line expression uses
bisect
to find the root.
>> r = bisect('myfx',xb)
r =
3.5214
We store the root in r
so that we can use that value later.
4. Verify that a Root Was Found
A good engineer always checks the result of a calculation. That check can be a simple "does-this-make-sense" test, or it can be an automated procedure for estimating the error in a calculation.
At the start of this root-finding example, we made a plot of .
From that plot we saw that crosses the axis between
and . Therefore, the value of 3.5214 returned by
bisect
seems like a good guess at the root.
Since the computational machinery is readily available, it always makes sense to verify the root returned by bisection. We simply evaluate
>> fr = myfx(r)
fr =
0
In this case fr
is exactly zero, which is surprising result. Normally
we expect the value of to be small, but not exactly zero.
We can create a visual confirmation of the root by adding the point
to the plot created by brackPlot
or to a plot created
without brackPlot
. Here we add the root to the plot created
by `brackPlot'.
Running these statements
>> hold('on') % Assumes the plot is still active
>> plot(r,fr,'o') % Add the guess at the root
>> hold('off')
produces a plot with the root indicated by a circle.
By plotting the suspected root, we confirm that is small compared to the value of the function. We also see that is the location where crosses the axis.
Complete MATLAB
Solution
A complete MATLAB
for this root-finding recipe is
is listed below
(download the code).
function fxroot
% fxroot Demonstrate basic root-finding recipe with brackPlot and bisect
% -- Get possible interval with the root
xb = brackPlot('myfx',0,5);
% -- Use bisection to refine the interval
r = bisect('myfx',xb);
% -- Add the value returned by bisect for visual confirmation
hold('on')
fr = myfx(r);
plot(r,fr,'*')
hold('off')
fprintf('\nRoot found at x = %8.5f where f(x) = %12.3e\n',r,fr)
end
Advice for using this recipe
-
First, write a simple m-file to evaluate the function
-
Plot to visually inspect its behavior. Not all are as well-behaved as the one used in this example.
-
Find a bracket to the root. You do not need to use
brackPlot
!brackPlot
is a convenience function that simplifies the search for brackets. It is not fool-proof. For root-finding problems embedded in more complex analysis, you will probably want to avoid usingbrackPlot
and instead develop a custom algorithm for guessing the roots. -
Always test your root to make sure
bisect
, or whatever root-finding algorithm you are using, has found a root, not a singularity.