2022-01-21 Rootfinding
Contents
2022-01-21 Rootfinding¶
Office hours: Monday 9-10pm, Tuesday 2-3pm, Thursday 2-3pm
Next week will stay virtual. Plan is to start in-person the following Monday (Jan 31)
Last time¶
More on floating point
Discuss Taylor Series activity
Condition numbers
Reliable = well-conditioned and stable¶
Mathematical functions \(f(x)\) can be ill-conditioned (big \(\kappa\))¶
Modeling is how we turn an abstract question into a mathematical function
We want well-conditioned models (small \(\kappa\))
Some systems are intrinsically sensitive: fracture, chaotic systems, combustion
Algorithms f(x) can be unstable¶
Unreliable, though sometimes practical
Unstable algorithms are constructed from ill-conditioned parts
An ill-conditioned problem from Paul Olum¶
From Surely You’re Joking, Mr. Feynman (page 113)
So Paul is walking past the lunch place and these guys are all excited. “Hey, Paul!” they call out. “Feynman’s terrific! We give him a problem that can be stated in ten seconds, and in a minute he gets the answer to 10 percent. Why don’t you give him one?” Without hardly stopping, he says, “The tangent of 10 to the 100th.” I was sunk: you have to divide by pi to 100 decimal places! It was hopeless.
What’s the condition number?
\(f(x) = \tan x\)
\(f'(x) = 1 + \tan^2 x\)
tan(1e100)
-0.4116229628832498
tan(BigFloat("1e100", precision=500))
0.4012319619908143541857543436532949583238702611292440683194415381168718098221194
Forward and backward error¶
plot(tan, ylims=(-10, 10), x=LinRange(-10, 10, 1000))
Read Backward error¶
For ill-conditioned functions, the best we can hope for is small backward error.
Feynman could have answered with any real number and it would have a (relative) backward error less than \(10^{-100}\). All the numbers we saw on the previous slide had tiny backward error because
for some \(\epsilon < 10^{-100}\). (All real numbers can be written this way for some small \(\epsilon\).)
Go find some functions…¶
Find a function \(f(x)\) that models something you’re interested in
Plot its condition number \(\kappa\) as a function of \(x\)
Plot the relative error (using single or double precision; compare using Julia’s
big)Is the relative error ever much bigger than \(\kappa \epsilon_{\text{machine}}\)?
Can you find what caused the instability?
Share on Zulip
Further reading: FNC Introduction¶
Stability demo: Volume of a polygon¶
X = [1 0; 2 1; 1 3; 0 1; -1 1.5; -2 -1; .5 -2; 1 0]
R(θ) = [cos(θ) -sin(θ); sin(θ) cos(θ)]
Y = X * R(deg2rad(0))' .+ [-1 0]
plot(Y[:,1], Y[:,2], seriestype=:shape, aspect_ratio=:equal, xlims=(-3, 3), ylims=(-3, 3))
using LinearAlgebra
function pvolume(X)
n = size(X, 1)
vol = sum(det(X[i:i+1, :] .- X[1:1, :]) / 2 for i in 1:n-1)
end
@show ref = pvolume(Y)
[det(Y[i:i+1, :]) for i in 1:size(Y, 1)-1]
ref = pvolume(Y) = 9.25
7-element Vector{Float64}:
0.0
3.0
3.0
0.5
6.5
5.5
-0.0
What happens if this polygon is translated, perhaps far away?
Rootfinding¶
Given \(f(x)\), find \(x\) such that \(f(x) = 0\).
We’ll work with scalars (\(f\) and \(x\) are just numbers) for now, and revisit later when they vector-valued.
Change inputs to outputs¶
\(f(x; b) = x^2 - b\)
\(x(b) = \sqrt{b}\)
\(f(x; b) = \tan x - b\)
\(x(b) = \arctan b\)
\(f(x) = \cos x + x - b\)
\(x(b) = ?\)
We aren’t given \(f(x)\), but rather an algorithm f(x) that approximates it.
Sometimes we get extra information, like
fp(x)that approximates \(f'(x)\)If we have source code for
f(x), maybe it can be transformed “automatically”
Example: Queueing¶
In a simple queueing model, there is an arrival rate and a departure (serving) rate. While waiting in the queue, there is a probability of “dropping out”. The length of the queue in this model is
One model for the waiting time (where these rates are taken from exponential distributions) is
wait(arrival, departure) = log(arrival / departure) / departure
plot(d -> wait( , d), xlims=(.1, 1), xlabel="departure", ylabel="wait")
syntax: unexpected ","
Stacktrace:
[1] top-level scope
@ In[7]:3
[2] eval
@ ./boot.jl:373 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1196
Departure rate given wait¶
Easy to measure wait
I have a limited tolerance for waiting
my_wait = 0.8
plot([d -> wait(1, d) - my_wait, d -> 0], xlims=(.1, 1), xlabel="departure", ylabel="wait")
using Roots
d0 = find_zero(d -> wait(1, d) - my_wait, 1)
scatter!([d0], [0], marker=(:circle), title="d0 = $d0")
Example: Nonlinear Elasticity¶
Strain-energy formulation¶
where \(I_1 = \lambda_1^2 + \lambda_2^2 + \lambda_3^3\) and \(J = \lambda_1 \lambda_2 \lambda_3\) are invariants defined in terms of the principle stretches \(\lambda_i\).
Uniaxial extension¶
In the experiment, we would like to know the stress as a function of the stretch \(\lambda_1\). We don’t know \(J\), and will have to determine it by solving an equation.
How much does the volume change?¶
Using symmetries of uniaxial extension, we can write an equation \(f(\lambda, J) = 0\) that must be satisfied. We’ll need to solve a rootfinding problem to compute \(J(\lambda)\).
function f(lambda, J)
C_1 = 1.5e6
D_1 = 1e8
D_1 * J^(8/3) - D_1 * J^(5/3) + C_1 / (3*lambda) * J - C_1 * lambda^2/3
end
plot(J -> f(4, J), xlims=(0.1, 3), xlabel="J", ylabel="f(J)")
find_J(lambda) = find_zero(J -> f(lambda, J), 1)
plot(find_J, xlims=(0.1, 5), xlabel="lambda", ylabel="J")
An algorithm: Bisection¶
Bisection is a rootfinding technique that starts with an interval \([a,b]\) containing a root and does not require derivatives. Suppose \(f\) is continuous. What is a sufficient condition for \(f\) to have a root on \([a,b]\)?
hasroot(f, a, b) = f(a) * f(b) < 0
f(x) = cos(x) - x
plot(f)
Bisection¶
function bisect(f, a, b, tol)
mid = (a + b) / 2
if abs(b - a) < tol
return mid
elseif hasroot(f, a, mid)
return bisect(f, a, mid, tol)
else
return bisect(f, mid, b, tol)
end
end
x0 = bisect(f, -1, 3, 1e-5)
x0, f(x0)
(0.7390861511230469, -1.7035832658995886e-6)
How fast does it converge?¶
function bisect_hist(f, a, b, tol)
mid = (a + b) / 2
if abs(b - a) < tol
return [mid]
elseif hasroot(f, a, mid)
return prepend!(bisect_hist(f, a, mid, tol), [mid])
else
return prepend!(bisect_hist(f, mid, b, tol), [mid])
end
end
bisect_hist (generic function with 1 method)
bisect_hist(f, -1, 3, 1e-4)
17-element Vector{Float64}:
1.0
0.0
0.5
0.75
0.625
0.6875
0.71875
0.734375
0.7421875
0.73828125
0.740234375
0.7392578125
0.73876953125
0.739013671875
0.7391357421875
0.73907470703125
0.739105224609375
Iterative bisection¶
Data structures often optimized for appending rather than prepending.
Bounds stack space
function bisect_iter(f, a, b, tol)
hist = Float64[]
while abs(b - a) > tol
mid = (a + b) / 2
push!(hist, mid)
if hasroot(f, a, mid)
b = mid
else
a = mid
end
end
hist
end
bisect_iter (generic function with 1 method)
bisect_iter(f, -1, 3, 1e-4)
16-element Vector{Float64}:
1.0
0.0
0.5
0.75
0.625
0.6875
0.71875
0.734375
0.7421875
0.73828125
0.740234375
0.7392578125
0.73876953125
0.739013671875
0.7391357421875
0.73907470703125
Let’s plot the error¶
where \(r\) is the true root, \(f(r) = 0\).
r = bisect(f, -1, 3, 1e-14) # what are we trusting?
hist = bisect_hist(f, -1, 3, 1e-10)
scatter( abs.(hist .- r), yscale=:log10)
ks = 1:length(hist)
plot!(ks, 4 * (.5 .^ ks))
Evidently the error \(e_k = x_k - x_*\) after \(k\) bisections satisfies the bound
Next time¶
Discuss rootfinding examples
Limitations of bisection
Convergence classes
Newton methods
convergence theory
derivation via convergence theory
reliability
For you¶
Share an equation that is useful for one task, but requires rootfinding for another task.
OR: explain a system in which one stakeholder has the natural inputs, but a different stakeholder only knows outputs.
Please collaborate.