Read about this topic here: Solving for zeros with julia.
For the impatient, these questions are related to the zeros of a real-valued function $f$. That is, those values of $x$ with $f(x)=0$.
Finding zeros of a polynomial (called "roots" when the function is a polynomial) is a familiar task that can be aided by a few key equations, such as the quadratic equation. However, in general, finding a zero of a function will require a numeric approach. The Roots
package of Julia
will provide some features. This is loaded when MTH229
is:
using MTH229
using Plots
plotly()
┌ Info: For saving to png with the Plotly backend PlotlyBase has to be installed. └ @ Plots /Users/verzani/.julia/packages/Plots/1KWPG/src/backends.jl:318
Plots.PlotlyBackend()
Graphically, a zero of a continuous function $f(x)$ occurs where the graph crosses or touches the $x$-axis. Without much work, a zero can be estimated to one or two decimal points from a graph. For example, we can zoom in on the zero of $f(x) = x^5 + x - 1$ by graphing over $[0,1]$:
f(x) = x^5 + x - 1
plot(f, 0, 1)
plot!(zero, 0, 1)
We can see that the zero is near $0.75$, but be careful reading too much into a graph. Since there are only so many pixels in a graph, and typically even fewer points chosen, what looks like a curve is really just a stick figure if you zoom in far enough. Replotting over a smaller domain can give more accuracy, but it is better to use a graph to get a sense of where the desired answer is and then use a numeric method to "zoom" in on the answer. In this project we discuss one such method for "zooming in" – the bisection method.
The notes mention the bisection method for zero-finding. This is based on the intermediate value theorem:
The intermediate value theorem states that if a continuous function, $f$, over an interval, $[a, b]$, takes values $f(a)$ and $f(b)$ at each end of the interval, then it also takes any value between $f(a)$ and $f(b)$ at some point within the interval.
For our purposes, this is specialized to "Bolzano's theorem":
If a continuous function has values of opposite sign inside a closed interval, then it has a zero in that interval
In particular, if $f(x)$ is continuous on $[a,b]$ and $f(a)$ and $f(b)$ have different signs, then there must be a value $c$ with $a < c < b$ where $f(c) = 0$.
There may be more than one zero, but there is a guarantee of at least one.
Not all functions can have their real zeros solved algebraically, and not all applications can be answered by the accuracy provided by a graph. In such situations, numeric methods may be of interest.
The MTH229
package defines a bisection
method implementing the bisection method, it assumes it has been passed values $a < b$ with $f(a)$ and $f(b)$ having different sign. In short, [a,b]
is a bracketing interval for $f$.
The algorithm to find a value $c$ proceeds in steps. At each step the bracketing interval is split in two at a value $m$. One – and only one – of these three things must be true: either $[a,m]$ is a bracketing interval, $[m,b]$ is a bracketing interval, or $f(m) = 0$. If the latter, the algorithm stops, otherwise the new bracketing interval splits the original one in half and the process proceeds again.
Running this demo illustrates.
bisection(sin, 3, 4)
a###############################################################b a###############################b................................ a###############b................................................ ........a#######b................................................ ........a###b.................................................... ........a#b...................................................... .........ab...................................................... ⋮⋮
3.141592653589793
In the demo, after $8$ iterations, there aren't enough pixels to see more subdivisions, but mathematically, unless the algorithm finds an exact zero, this process would continue infinitely, with the bracketing interval getting infinitely small. In the process this traps the zero.
On the computer, the process basically stops when the size of the bracketing interval gets too small to subdivide using floating point numbers, unless instructed otherwise.
In the Roots
package is the fzero
method that implements the bisection method, only a bit more carefully. The MTH229
package loads this for you.
For a bracketing interval, it is guaranteed to find a c
such that the function changes sign between adjacent floating point values around c
, or c
is an exact zero. It is used as: fzero(f, a, b)
:
f(x) = x^2 - 2
fzero(f, 1, 2) # finds sqrt(2)
1.4142135623730951
Many problems are more naturally expressed by solving $f(x) = g(x)$, and not $f(x) = 0$, as expected by fzero
. This is no issue, as it only requires the extra step of defining a function of the differences $h(x) = f(x) - g(x)$, as $h(x) = 0$ implies $f(x) = g(x)$.
For example, consider this question:
Find the intersection point of $4 - e^{x/10} = e^{x/15}$ by first graphing to see approximately where the answer is. From the graph, identify a bracket and then use fzero
to numerically estimate the intersection point.
We could plot both functions:
f(x) = 4 - exp(x/10)
g(x) = exp(x/15)
plot(f, 0, 20)
plot!(g, 0, 20)
Or we could plot the difference:
h(x) = f(x) - g(x)
plot(h, 0, 20)
plot!(zero, 0, 20)
From either graph, we see quickly that the interval $[5,10]$ will be a bracketing interval for $h$, so we can find the intersection point with:
fzero(h, 5, 10)
8.205886667065423
How to find all zeros of a function within a given interval? The bisection method guarantees only one between a bracketing interval. Well, a simple algorithm–-which is not guaranteed to find all zeros, but should do pretty well–-is to divide the interval into many subintervals and for each, see if the subinterval is a bracketing interval, and if so find a zero in that subinterval.
The fzeros(f, a, b)
function does basically that, though it works a bit harder to find answers even when the subinterval is not a bracketing interval. It makes no guarantees though, as the bisection method can–-and may–-miss some answers. The fzeros
function is only pretty good at finding all the zeros over the interval. For some functions – especially those cooked up by clever math professors – the choice of interval can lead to fzeros
finding fewer than is mathematically possible. The function should be used in combination with a plot and with as narrow an interval specified, as reasonable.
Since fzeros
can return 0, 1, or more zeros, it uses a container to hold its answers. This means that some numbers won't print with 16 digits displayed. This does not mean the accuracy is lost, it is just not displayed.
The output of fzeros
is a collection of values. It may be desirable to pass these onto another function. This is essentially composition. For example, we can check our work using this pattern:
f(x) = cos(x) + cos(2x)
zs = fzeros(f, 0, 2pi)
f.(zs) # broadcast f over each of the zs
3-element Vector{Float64}: 3.3306690738754696e-16 0.0 -1.9984014443252818e-15
We see that the values are all basically $0$, save for round-off error.
The key above is the use of the "dot" before the (
to broadcast f
over each value in the container zs
. Just entering f(zs)
would lead to error about no matching method.