Chapter 4

using Plots, SpecialFunctions
J₃(x) = besselj(3, x)
plot(J₃, 0, 20;
    title="Bessel function",
    xaxis=(L"x"),  yaxis=(L"J_3(x)"),  grid=:xy)

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.

using NLsolve
ω = []
for guess = [6., 10. ,13., 16., 19.]
    s = nlsolve(x -> J₃(x[1]), [guess], ftol=1e-14)
    append!(ω, s.zero)
end

y = J₃.(ω)
@pt :header=["root estimate", "function value"] [ω y]

scatter!(ω, y, 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$.

r = []
for guess = [3., 6., 10., 13.]
    f(x) = J₃(x[1]) - 0.2
    s = nlsolve(f, [guess], ftol=1e-14)
    append!(r, s.zero)
end
scatter!(r, J₃.(r), title="Roots and other Bessel values")

Consider first the function

f(x) = (x - 1) * (x - 2);

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.

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

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:

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!

Let’s convert the roots of a quadratic polynomial $f(x)$ to a fixed point problem.

using Polynomials
p = Polynomial([3.5, -4,1])
r = roots(p)
rmin, rmax = extrema(r)
@show rmin, rmax;

(rmin, rmax) = (1.2928932188134525, 2.7071067811865475)

We define $g(x)=x-p(x)$.

g(x) = x - p(x)

Intersections of $y=g(x)$ with the line $y=x$ are fixed points of $g$ and thus roots of $f$. (Only one is shown in the chosen plot range.)

using Plots
plt = plot([g x->x], 2, 3;
    l=2, label=[L"y=g(x)" L"y=x"],
    xlabel=L"x",  ylabel=L"y", 
    aspect_ratio=1,
    title="Finding a fixed point",  legend=:bottomright)

If we evaluate $g(2.1)$, we get a value of almost 2.6, so this is not a fixed point.

x = 2.1;
y = g(x)

However, $y=g(x)$ is considerably closer to the fixed point at around 2.7 than $x$ is. Suppose then that we adopt $y$ as our new $x$ value. Changing the $x$ coordinate in this way is the same as following a horizontal line over to the graph of $y=x$.

plot!([x, y], [y, y], arrow=true, color=3)

Now we can compute a new value for $y$. We leave $x$ alone here, so we travel along a vertical line to the graph of $g$.

x = y;  y = g(x)
plot!([x, x], [x, y], arrow=true, color=4)

You see that we are in a position to repeat these steps as often as we like. Let’s apply them a few times and see the result.

for k = 1:5
    plot!([x, y], [y, y], color=3);  
    x = y       # y becomes the new x
    y = g(x)    # g(x) becomes the new y
    plot!([x, x], [x, y], color=4)  
end
plt

The process spirals in beautifully toward the fixed point we seek. Our last estimate has almost 4 accurate digits.

abs(y - rmax) / rmax

Now let’s try to find the other fixed point $\approx 1.29$ in the same way. We’ll use 1.3 as a starting approximation.

plt = plot([g x->x], 1, 2, l=2, label=["y=g(x)" "y=x"], aspect_ratio=1, 
    xlabel=L"x", ylabel=L"y", title="Divergence", legend=:bottomright)

x = 1.3; y = g(x);
arrow = false
for k = 1:5
    plot!([x, y], [y, y], arrow=arrow, color=3)  
    x = y       # y --> new x
    y = g(x)    # g(x) --> new y
    plot!([x, x], [x, y], arrow=arrow, color=4)
    if k > 2; arrow = true; end
end
plt

This time, the iteration is pushing us away from the correct answer.

We revisit Demo 4.2.1 and investigate the observed convergence more closely. Recall that above we calculated $g'(p)\approx-0.42$ at the convergent fixed point.

p = Polynomial([3.5, -4, 1])
r = roots(p)
rmin, rmax = extrema(r)
@show rmin, rmax;

(rmin, rmax) = (1.2928932188134525, 2.7071067811865475)

Here is the fixed-point iteration. This time we keep track of the whole sequence of approximations.

g(x) = x - p(x)
x = [2.1]
for k = 1:12
    push!(x, g(x[k]))
end
x

It’s illuminating to construct and plot the sequence of errors.

err = @. abs(x - rmax)
plot(0:12, err;
    m=:o,
    xaxis=("iteration number"),  yaxis=("error", :log10),
    title="Convergence of fixed-point iteration")

It’s quite clear that the convergence quickly settles into a linear rate. We could estimate this rate by doing a least-squares fit to a straight line. Keep in mind that the values for small $k$ should be left out of the computation, as they don’t represent the linear trend.

y = log.(err[5:12])
p = Polynomials.fit(5:12, y, 1)

We can exponentiate the slope to get the convergence constant σ.

σ = exp(p.coeffs[2])

The error should therefore decrease by a factor of σ at each iteration. We can check this easily from the observed data.

[err[i+1] / err[i] for i in 8:11]

The methods for finding σ agree well.

Suppose we want to find a root of the function

f(x) = x * exp(x) - 2
using Plots
plot(f, 0, 1.5; 
    label="function",  legend=:topleft,
    grid=:y,  ylim=[-2, 4],  xlabel=L"x",  ylabel=L"y")

From the graph, it is clear that there is a root near $x=1$. So we call that our initial guess, $x_1$.

x₁ = 1
y₁ = f(x₁)
scatter!([x₁], [y₁], label="initial point")

Next, we can compute the tangent line at the point $\bigl(x_1,f(x_1)\bigr)$, using the derivative.

df_dx(x) = exp(x) * (x + 1)
m₁ = df_dx(x₁)
tangent = x -> y₁ + m₁ * (x - x₁)

plot!(tangent, 0, 1.5, l=:dash, label="tangent line",
    title="Tangent line approximation")

In lieu of finding the root of $f$ itself, we settle for finding the root of the tangent line approximation, which is trivial. Call this $x_2$, our next approximation to the root.

@show x₂ = x₁ - y₁ / m₁
scatter!([x₂], [0], label="tangent root", title="First iteration")

y₂ = f(x₂)

The residual (i.e., value of $f$) is smaller than before, but not zero. So we repeat the process with a new tangent line based on the latest point on the curve.

plot(f, 0.82, 0.87;
    label="function",  legend=:topleft,
    xlabel=L"x",  ylabel=L"y",
    title="Second iteration")

scatter!([x₂], [y₂], label="starting point")

m₂ = df_dx(x₂)
tangent = x -> y₂ + m₂ * (x - x₂)
plot!(tangent, 0.82, 0.87; l=:dash, label="tangent line")

@show x₃ = x₂ - y₂ / m₂
scatter!([x₃], [0], label="tangent root")

y₃ = f(x₃)

Judging by the residual, we appear to be getting closer to the true root each time.

We again look at finding a solution of $x e^x=2$ near $x=1$. To apply Newton’s method, we need to calculate values of both the residual function $f$ and its derivative.

f(x) = x * exp(x) - 2;
df_dx(x) = exp(x) * (x + 1);

We don’t know the exact root, so we use nlsolve to determine a proxy for it.

using NLsolve
r = nlsolve(x -> f(x[1]), [1.0]).zero

We use $x_1=1$ as a starting guess and apply the iteration in a loop, storing the sequence of iterates in a vector.

x = [1; zeros(4)]
for k = 1:4
    x[k+1] = x[k] - f(x[k]) / df_dx(x[k])
end
x

Here is the sequence of errors.

ϵ = @. x - r

Because the error reaches machine epsilon so rapidly, we’re going to use extended precision to allow us to take a few more iterations. We’ll take the last iteration as the most accurate root estimate.

x = [BigFloat(1); zeros(7)]
for k = 1:7
    x[k+1] = x[k] - f(x[k]) / df_dx(x[k])
end
r = x[end]

ϵ = @. Float64(x[1:end-1] - r)

The exponents in the scientific notation definitely suggest a squaring sequence. We can check the evolution of the ratio in (4.3.9).

logerr = @. log10(abs(ϵ))
ratios = [NaN; [logerr[i+1] / logerr[i] for i in 1:length(logerr)-1]]
@pt :header=["iteration", "error", "log error", "ratio"] [1:7 ϵ logerr ratios]

The clear convergence to 2 above constitutes good evidence of quadratic convergence.

Suppose we want to evaluate the inverse of the function $h(x)=e^x-x$. This means solving $y=e^x-x$ for $x$ when $y$ is given, which has no elementary form. If a value of $y$ is given numerically, though, we simply have a rootfinding problem for $f(x)=e^x-x-y$.

g(x) = exp(x) - x
dg_dx(x) = exp(x) - 1
y = range(g(0), g(2), 200)
x = zeros(length(y))
for (i, y) in enumerate(y)
    f(x) = g(x) - y
    df_dx(x) = dg_dx(x)
    r = FNC.newton(f, df_dx, y)
    x[i] = r[end]
end

plot(g, 0, 2, aspect_ratio=1, label=L"g(x)")
plot!(y, x, label=L"g^{-1}(y)", title="Function and its inverse")
plot!(x -> x, 0, maximum(y), label="", l=(:dash, 1), color=:black)

We return to finding a root of the equation $x e^x=2$.

using Plots
f(x) = x * exp(x) - 2;
plot(f, 0.25, 1.25;
    label="function",  legend=:topleft,
    xlabel=L"x",  ylabel=L"y")

From the graph, it’s clear that there is a root near $x=1$. To be more precise, there is a root in the interval $[0.5,1]$. So let us take the endpoints of that interval as two initial approximations.

x₁ = 1;
y₁ = f(x₁);
x₂ = 0.5;
y₂ = f(x₂);
scatter!([x₁, x₂], [y₁, y₂];
    label="initial points",
    title="Two initial values")

Instead of constructing the tangent line by evaluating the derivative, we can construct a linear model function by drawing the line between the two points $\bigl(x_1,f(x_1)\bigr)$ and $\bigl(x_2,f(x_2)\bigr)$. This is called a secant line.

m₂ = (y₂ - y₁) / (x₂ - x₁)
secant = x -> y₂ + m₂ * (x - x₂)
plot!(secant, 0.25, 1.25, label="secant line", l=:dash, color=:black,
    title="Secant line")

As before, the next root estimate in the iteration is the root of this linear model.

x₃ = x₂ - y₂ / m₂
@show y₃ = f(x₃)
scatter!([x₃], [0], label="root of secant", title="First iteration")

For the next linear model, we use the line through the two most recent points. The next iterate is the root of that secant line, and so on.

m₃ = (y₃ - y₂) / (x₃ - x₂)
x₄ = x₃ - y₃ / m₃

We check the convergence of the secant method from Demo 4.4.1. Again we will use extended precision to get a longer sequence than double precision allows.

f(x) = x * exp(x) - 2
x = FNC.secant(f, BigFloat(1), BigFloat(0.5), xtol=1e-80, ftol=1e-80);

We don’t know the exact root, so we use the last value as a proxy.

r = x[end]

Here is the sequence of errors.

ϵ = @. Float64(r - x[1:end-2])

It’s not easy to see the convergence rate by staring at these numbers. We can use (4.4.8) to try to expose the superlinear convergence rate.

logerr = @. log10(abs(ϵ))
ratios = [NaN; [logerr[i+1] / logerr[i] for i in 1:length(logerr)-1]]
@pt :header=["iteration", "error", "log error", "ratio"] [eachindex(ϵ) ϵ logerr ratios]

As expected, this settles in at around 1.618.

Here we look for a root of $x+\cos(10x)$ that is close to 1.

f(x) = x + cos(10 * x)
interval = [0.5, 1.5]

plot(f, interval..., label="Function", legend=:bottomright,
    grid=:y, ylim=[-0.1, 3], xlabel=L"x", ylabel=L"y")

We choose three values to get the iteration started.

x = [0.8, 1.2, 1]
y = @. f(x)
scatter!(x, y, label="initial points")

If we were using forward interpolation, we would ask for the polynomial interpolant of $y$ as a function of $x$. But that parabola has no real roots.

using Polynomials
q = Polynomials.fit(x, y, 2)      # interpolating polynomial
plot!(x -> q(x), interval..., l=:dash, label="interpolant")

To do inverse interpolation, we swap the roles of $x$ and $y$ in the interpolation.

```{tip}
:class: dropdown
By giving two functions in the plot call, we get the parametric plot $(q(y),y)$ as a function of $y$.
```

```{code-cell}
plot(f, interval..., label="Function",
    legend=:bottomright, grid=:y, xlabel=L"x", ylabel=L"y")
scatter!(x, y, label="initial points")

q = Polynomials.fit(y, x, 2)       # interpolating polynomial
plot!(y -> q(y), y -> y, -0.1, 2.6, l=:dash, label="inverse interpolant")
```

We seek the value of $x$ that makes $y$ zero. This means evaluating $q$ at zero.

```{code-cell}
q(0)
```

Let's restart the process with `BigFloat` numbers to get a convergent sequence.

```{code-cell}
x = BigFloat.([8, 12, 10]) / 10
y = @. f(x)

for k = 3:12
    q = Polynomials.fit(y[k-2:k], x[k-2:k], 2)
    push!(x, q(0))
    push!(y, f(x[k+1]))
end

println("residual = $(f(x[end]))")
```

As far as our current precision is concerned, we have an exact root.

```{code-cell}
r = x[end]
ϵ = @. Float64(abs(r - x[1:end-1]))
logerr = @. log10(abs(ϵ))
ratios = [NaN; [logerr[i+1] / logerr[i] for i in 1:length(logerr)-1]]
@pt :header=["iteration", "error", "log error", "ratio"] [eachindex(ϵ) ϵ logerr ratios]
```

The convergence is probably superlinear at a rate of $\alpha=1.8$ or so.

A system of nonlinear equations is defined by its residual and Jacobian.

function func(x)
    [exp(x[2] - x[1]) - 2,
        x[1] * x[2] + x[3],
        x[2] * x[3] + x[1]^2 - x[2]
    ]
end;

function jac(x)
    [
        -exp(x[2] - x[1]) exp(x[2] - x[1]) 0
        x[2] x[1] 1
        2*x[1] x[3]-1 x[2]
    ]
end;

We will use a BigFloat starting value, and commensurately small stopping tolerances, in order to get a sequence long enough to measure convergence.

x₁ = BigFloat.([0, 0, 0])
ϵ = eps(BigFloat)
x = FNC.newtonsys(func, jac, x₁, xtol=ϵ, ftol=ϵ);

Let’s compute the residual of the last result in order to check the quality.

r = x[end]
@show residual = norm(func(r));

residual = norm(func(r)) = 0.0

We take the sequence of norms of errors, applying the log so that we can look at the exponents.

logerr = [Float64(log(norm(r - x[k]))) for k in 1:length(x)-1]
ratios = [NaN; [logerr[i+1] / logerr[i] for i in 1:length(logerr)-1]]
@pt :header=["iteration", "log error", "ratio"] [eachindex(logerr) logerr ratios]

The ratio is neatly converging toward 2, which is expected for quadratic convergence.

To solve a nonlinear system, we need to code only the function defining the system, and not its Jacobian.

f(x) = 
    [
        exp(x[2] - x[1]) - 2,
        x[1] * x[2] + x[3],
        x[2] * x[3] + x[1]^2 - x[2]
    ]

In all other respects usage is the same as for the newtonsys function.

x₁ = [0.0, 0.0, 0.0]
x = FNC.levenberg(f, x₁)

It’s always a good idea to check the accuracy of the root, by measuring the residual (backward error).

r = x[end]
println("backward error = $(norm(f(r)))")

backward error = 1.2707848769787674e-13

Looking at the convergence in norm, we find a convergence rate between linear and quadratic, like with the secant method.

logerr = [log(norm(r - x[k])) for k in 1:length(x)-1]
ratios = [NaN; [logerr[i+1] / logerr[i] for i in 1:length(logerr)-1]]
@pt :header=["iteration", "log error", "ratio"] [eachindex(logerr) logerr ratios]

We will observe the convergence of Function 4.6.3 for different levels of the minimum least-squares residual. We start with a function mapping from $\real^2$ into $\real^3$, and a point that will be near the optimum.

g(x) = [sin(x[1] + x[2]), cos(x[1] - x[2]), exp(x[1] - x[2])]
p = [1, 1];

The function $\mathbf{g}(\mathbf{x}) - \mathbf{g}(\mathbf{p})$ obviously has a zero residual at $\mathbf{p}$. We’ll make different perturbations of that function in order to create nonzero residuals.

using Printf
plt = plot(xlabel="iteration", yaxis=(:log10, "error"),
    title="Convergence of Gauss–Newton")
for R in [1e-3, 1e-2, 1e-1]
    # Define the perturbed function.
    f(x) = g(x) - g(p) + R * normalize([-1, 1, -1])
    x = FNC.levenberg(f, [0, 0])
    r = x[end]
    err = [norm(x - r) for x in x[1:end-1]]
    normres = norm(f(r))
    plot!(err, label=@sprintf("R=%.2g", normres))
end
plt

In the least perturbed case, where the minimized residual is less than 10^-3, the convergence is plausibly quadratic. At the next level up, the convergence starts similarly but suddenly stagnates for a long time. In the most perturbed case, the quadratic phase is nearly gone and the overall shape looks linear.

m = 25;
s = range(0.05, 6, length=m)
ŵ = @. 2 * s / (0.5 + s)                      # exactly on the curve
w = @. ŵ + 0.15 * cos(2 * exp(s / 16) * s);     # smooth noise added

scatter(s, w, label="noisy data",
    xlabel="s", ylabel="v", leg=:bottomright)
plot!(s, ŵ, l=:dash, color=:black, label="perfect data")

The idea is to pretend that we know nothing of the origins of this data and use nonlinear least squares to find the parameters in the theoretical model function $v(s)$. In (4.7.2), the $s$ variable plays the role of $t$, and $v$ plays the role of $g$.

function misfit(x)
    V, Km = x   # rename components for clarity
    return @. V * s / (Km + s) - w
end

In the Jacobian the derivatives are with respect to the parameters in $\mathbf{x}$.

function misfitjac(x)
    V, Km = x   # rename components for clarity
    J = zeros(m, 2)
    J[:, 1] = @. s / (Km + s)              # dw/dV
    J[:, 2] = @. -V * s / (Km + s)^2         # dw/d_Km
    return J
end

x₁ = [1, 0.75]
x = FNC.newtonsys(misfit, misfitjac, x₁)

@show V, Km = x[end]  # final values

(V, Km) = x[end] = [1.968652598378229, 0.4693037307416775]

The final values are reasonably close to the values $V=2$, $K_m=0.5$ that we used to generate the noise-free data. Graphically, the model looks close to the original data.

model(s) = V * s / (Km + s)
plot!(model, 0, 6, label="nonlinear fit")

For this particular model, we also have the option of linearizing the fit process. Rewrite the model as

for the new fitting parameters $\alpha=K_m/V$ and $\beta=1/V$. This corresponds to the misfit function whose entries are

for $i=1,\ldots,m$. Although this misfit is nonlinear in $s$ and $w$, it’s linear in the unknown parameters α and β. This lets us pose and solve it as a linear least-squares problem.

A = [s .^ (-1) s .^ 0]
u = 1 ./ w
α, β = A \ u

The two fits are different because they do not optimize the same quantities.

linmodel(x) = 1 / (β + α / x)
plot!(linmodel, 0, 6, label="linearized fit")

The truly nonlinear fit is clearly better in this case. It optimizes a residual for the original measured quantity rather than a transformed one we picked for algorithmic convenience.