```
```{raw} latex
%%start demo%%
```
In the theory of vibrations of a circular drum, the displacement of the drumhead can be expressed in terms of pure harmonic modes,
$$J_m(\omega_{k,m} r) \cos(m\theta) \cos(c \omega_{k,m} t),$$
where $(r,\theta)$ are polar coordinates, $0\le r\le 1$, $t$ is time, $m$ is a positive integer, $c$ is a material parameter, and $J_m$ is a _Bessel function of the first kind_. The quantity $\omega_{k,m}$ is a resonant frequency and is a positive root of the equation
$$J_m(\omega_{k,m}) = 0,$$
which states that the drumhead is clamped around the rim. Bessel functions often appear in physical problems featuring radial symmetry, and tabulating approximations to the zeros of Bessel functions occupied numerous mathematician-hours before computers were on the scene.
```{code-cell}
J3(x) = besselj(3,x)
plot(J3,0,20,title="Bessel function",
xaxis=(L"x"),yaxis=(L"J_3(x)"),grid=:xy)
```
::::{panels}
:column: col-7 left-side
:card: border-0 shadow-none
```{raw} latex
\begin{minipage}[t]{0.5\textwidth}
```
From the graph we see roots near 6, 10, 13, 16, and 19. We use `nlsolve` from the `NLsolve` package to find these roots accurately. It uses vector variables, so we have to code accordingly.
```{raw} latex
\end{minipage}\hfill
```
---
:column: col-5 right-side
:card: shadow-none comment
```{raw} latex
\begin{minipage}[t]{0.4\textwidth}\begin{mdframed}[default]\small
```
Type `\omega` followed by `Tab` to get the character `ω`.
The argument `ftol=1e-14` below is called a **keyword argument**. Here it sets a goal for the maximum value of $|f(x)|$.
```{raw} latex
\end{mdframed}\end{minipage}
```
::::
```{code-cell}
ω = []
for guess = [6.,10.,13.,16.,19.]
s = nlsolve(x->besselj(3,x[1]),[guess],ftol=1e-14)
append!(ω,s.zero)
end
```
```{code-cell}
pretty_table([ω J3.(ω)],["root estimate" "function value"])
```
```{code-cell}
scatter!(ω,J3.(ω),title="Bessel function with roots")
```
If instead we seek values at which $J_3(x)=0.2$, then we must find roots of the function $J_3(x)-0.2$.
```{code-cell}
r = []
for guess = [3.,6.,10.,13.]
f = x -> besselj(3,x[1])-0.2
s = nlsolve(f,[guess],ftol=1e-14)
append!(r,s.zero)
end
scatter!(r,J3.(r),title="Roots and other Bessel values")
```
```{raw} html

```
```{raw} latex
%%end demo%%
```
## Conditioning, error, and residual
In the rootfinding problem, the data is a continuous function $f$ and the result is a root. (This overrides our Chapter 1 notation of $f$ as the map from data to result.) How does the result change in response to perturbations in $f$? We will compute an absolute condition number rather than a relative one.
You might wonder about the relevance of perturbing a function as data to a problem. If nothing else, the values of $f$ will be represented in floating point and thus subject to rounding error. Furthermore, in many applications, $f$ might not be a simple formula but the result of a computation that uses an inexact algorithm. While there are infinitely many possible perturbations to a function, a constant perturbation is enough to get the main idea.
We assume $f$ has at least one continuous derivative near a particular root $r$. Suppose that $f$ is perturbed to $\tilde{f}(x) = f(x) + \epsilon$. As a result, the root (if it still exists) will be perturbed to $\tilde{r} = r + \delta$ such that $\tilde{f}(\tilde{r})=0$. We now compute an absolute condition number $\kappa_r$, which is the ratio $\left | \frac{\delta}{\epsilon} \right|$ as $\epsilon\to 0$.
Using Taylor's theorem,
```{math}
0 = f(r+\delta) + \epsilon \approx f(r) + f'(r) \delta + \epsilon.
```
Since $r$ is a root, we have $f(r)=0$. This lets us relate $\delta$ to $\epsilon$, and their ratio is the condition number.
```{index} condition number; of rootfinding
```
````{proof:theorem} Condition number of rootfinding
If $f$ is differentiable at a root $r$, then the absolute condition number of $r$ with respect to constant changes in $f$ is
```{math}
:label: rootcondnum
\kappa_r = \bigl| f'(r) \bigr|^{-1}.
```
We say $\kappa_r = \infty$ if $f'(r)=0$.
````
Equivalently, {eq}`rootcondnum` is just the magnitude of the derivative of the inverse function $f^{-1}$ at zero.
(demo-roots-cond)=
```{proof:demo}
```
```{raw} html
```
```{raw} latex
%%start demo%%
```
Consider first the function
```{code-cell}
f = x -> (x-1)*(x-2);
```
:::{index} ! Julia; splatting
:::
::::{panels}
:column: col-7 left-side
:card: border-0 shadow-none
```{raw} latex
\begin{minipage}[t]{0.5\textwidth}
```
At the root $r=1$, we have $f'(r)=-1$. If the values of $f$ were perturbed at every point by a small amount of noise, we can imagine finding the root of the function drawn with a thick ribbon, giving a range of potential roots.
```{raw} latex
\end{minipage}\hfill
```
---
:column: col-5 right-side
:card: shadow-none comment
```{raw} latex
\begin{minipage}[t]{0.4\textwidth}\begin{mdframed}[default]\small
```
The syntax `interval...` is called **splatting** and means to insert all the individual elements of `interval` as a sequence.
```{raw} latex
\end{mdframed}\end{minipage}
```
::::
```{code-cell}
interval = [0.8,1.2]
plot(f,interval...,ribbon=0.03,aspect_ratio=1,
xlabel=L"x",yaxis=(L"f(x)",[-0.2,0.2]))
scatter!([1],[0],title="Well-conditioned root")
```
The possible values for a perturbed root all lie within the interval where the ribbon intersects the $x$-axis. The width of that zone is about the same as the vertical thickness of the ribbon.
By contrast, consider the function
```{code-cell}
f = x -> (x-1)*(x-1.01);
```
Now $f'(1)=-0.01$, and the graph of $f$ will be much shallower near $x=1$. Look at the effect this has on our thick rendering:
```{code-cell}
plot(f,interval...,ribbon=0.03,aspect_ratio=1,
xlabel=L"x",yaxis=(L"f(x)",[-0.2,0.2]))
scatter!([1],[0],title="Poorly-conditioned root")
```
The vertical displacements in this picture are exactly the same as before. But the potential _horizontal_ displacement of the root is much wider. In fact, if we perturb the function entirely upward by the amount drawn here, the root disappears!
```{raw} html

```
```{raw} latex
%%end demo%%
```
```{index} ! residual; of rootfinding
```
We must accept that when $|f'|$ is small at the root, it may not be possible to get a small error in a computed root estimate. As always, the error is not a quantity we can compute without knowing the exact answer. There is something else we can measure, though.
::::{proof:definition} Rootfinding residual
If $\tilde{r}$ approximates a root $r$ of function $f$, then the **residual** at $\tilde{r}$ is $f(\tilde{r})$.
::::
It stands to reason that a small residual might be associated with a small error. To quantify the relationship, let $\tilde{r}$ approximate root $r$, and define the new function $g(x)=f(x)-f(\tilde{r})$. Trivially, $g(\tilde{r})=0$, meaning that $\tilde{r}$ is a true root of $g$. Since the difference $g(x)-f(x)$ is the residual value $f(\tilde{r})$, the residual is the distance to an exactly solved rootfinding problem.
```{index} backward error
```
```{proof:observation}
The backward error in a root estimate is equal to the residual.
```
In general, it is not realistic to expect a small error in a root approximation if the condition number {eq}`rootcondnum` is large. However, we can gauge the backward error from the residual.
## Multiple roots
```{index} ! roots; multiplicity of
```
The condition number {eq}`rootcondnum` naturally leads to the question of what happens if $f'(r)=0$ at a root $r$. The following definition agrees with and extends the notion of algebraic multiplicity in a polynomial to roots of more general differentiable functions.
(definition-rootproblem-simple)=
::::{proof:definition} Multiplicity of a root
If $f(r)=f'(r)=\cdots=f^{(m-1)}(r)=0$, but $f^{(m)}(r)\neq 0$, then we say $f$ has a root of **multiplicity** $m$ at $r$. In particular, if $f(r)=0$ and $f'(r)\neq 0$, then $m=1$ and we call $r$ a **simple root**.
::::
Another useful characterization of multiplicity $m$ is that $f(x)=(x-r)^m q(x)$ for a differentiable $q$ with $q(r)\neq 0$.
When $r$ is a nonsimple root, the condition number {eq}`rootcondnum` is effectively infinite.[^infcond] However, even if $r$ is simple, we should expect difficulty in rootfinding if the condition number is very large. This occurs when $|f'(r)|$ is very small, which means that the quotient $q$ satisfies $q(r)\approx 0$ and another root of $f$ is very close to $r$. We made the same observation about polynomial roots all the way back in {numref}`Demo {number}