<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Jupyter notebooks
This is a Jupyter notebook using Python.  You can install Jupyter locally to edit and interact with this notebook.
Rootfinding
Rootfinding is the process of solving $$f(x) = 0$$ for $x$.  The standard assumption is that $f 
Step1:  Which of these functions have at least one root?
Which have more than one root?
Can we determine these properties merely by evaluating $f(x)$ for some values of $x$?
Bisection
Bisection is a rootfinding technique that starts with an interval $[a,b]$ containing a root and does not require derivatives.
Step2:  Notice that we need to define hasroot above.
Let's try running it
Step3:  We get about 5 digits of accuracy.  Why?  How fast did we get there?
Step4:  Can you find any problems with this implementation?  List them below
Step5:  What's going wrong here?  How can we improve the implementation and what are fundamental limitations of the algorithm?
Convergence rate
Let's quantitatively revisit the convergence rate.  A convergent rootfinding algorithm produces a sequence of approximations $x_i$ such that $$\lim_{i \to \infty} x_i \to x_$$ where $f(x_) = 0$.  For analysis, it is convenient to define the errors $e_i = x_i - x_$. We say that an algorithm is linearly convergent* if $$\lim_{i \to \infty} |e_{i+1}| / |e_i| = \rho < 1.$$  A smaller convergence factor $\rho$ represents faster convergence.
What is $\rho$ for bisection?
Remarks on bisection
Specifying an interval is often inconvenient
An interval in which the function changes sign guarantees convergence (robustness)
No derivative information is required
Roots of even degree are problematic
The solution error is directly available
The convergence rate is modest -- one iteration per bit of accuracy
Newton-Raphson Method
Much of numerical analysis reduces to Taylor series, the approximation
$$ f(x) = f(x_0) + f'(x_0) (x-x_0) + f''(x_0) (x - x_0)^2 / 2 + \dotsb $$
centered on some reference point $x_0$.
In numerical computation, it is exceedingly rare to look beyond the first-order approximation
$$ \tilde f_{x_0}(x) = f(x_0) + f'(x_0)(x - x_0) . $$
Since $\tilde f_{x_0}(x)$ is a linear function, we can explicitly compute the unique solution of $\tilde f_{x_0}(x) = 0$ as
$$ x = x_0 - f(x_0) / f'(x_0) . $$
This is Newton's Method (aka Newton-Raphson or Newton-Raphson-Simpson) for finding the roots of differentiable functions.
Step6:  Oops, how can we fix this?
This output is kinda hard to read, so let's make it cleaner.
Step7:  Did we solve all of these equations?
How can the iteration break down?
Does choosing a different initial guess lead to different solutions?
How is this convergence test different from the one we used for bisection?
Is the convergence rate similar for all test equations?
Convergence of Newton-type algorithms
We would like to know sharp conditions on when Newton-type algorithms converge, and if so, how fast.  This theory will build on that for a general Fixed Point Iteration $x_{i+1} = g(x_i)$ where $g$ is a continuously differentiable function.  Suppose that there exists a fixed point $r = g(r)$.  By the mean value theorem, we have that
$$ x_{i+1} - r = g(x_i) - g(r) = g'(c_i) (x_i - r) $$
for some $c$ with $|c - r| < |x_i - r|$.
In other words, $|e_{i+1}| = |g'(c_i)| |e_i|$, which converges to zero if $|g'(c_i)| < 1$.
If $|g'(r)| < 1$ then for any $\epsilon > 0$ there is a neighborhood of $r$ such that $|g'(c)| < |g'(r)| + \epsilon$ for all $c$ in that neighborhood.
Consequently, we have
Step8:  We only get three digits correct despite a very small residual (and it takes many iterations to get there).
Difficulty computing zeros of polynomials can also arise when all the roots are simple.  For example, the Wilkinson polynomial
$$ \prod_{i=1}^{20} (x - i) = \sum_{i=0}^{20} a_i x^i $$
has roots at each of the positive integers up to 20, but the roots are extremely sensitive to perturbations of the coefficients $a_i$, as shown in this figure from Trefethen and Bau (1999).
Numerical difficulties in which "correct" algorithms produce unreliable solutions almost always stem from lack of stability and/or ill conditioning.
Absolute condition number
Consider a function $f
Step9:  Relative condition number
Given the relative nature of floating point arithmetic, it is more useful to discuss relative condition number,
$$ \kappa = \max_{\delta x} \frac{|\delta f|/|f|}{|\delta x|/|x|}
= \max_{\delta x} \Big[ \frac{|\delta f|/|\delta x|}{|f| / |x|} \Big] $$
or, if $f$ is differentiable,
$$ \kappa = \max_{\delta x} |f'(x)| \frac{|x|}{|f|} . $$
How does a condition number get big?
Take-home message
The relative accuracy of the best-case algorithm will not be reliably better than $\epsilon_{\text{machine}}$ times the condition number.
$$ \max_{\delta x} \frac{|\delta f|}{|f|} \ge \kappa \cdot \epsilon_{\text{machine}} $$
Numerical differentiation
Suppose we want to apply Newton's method to a function that we know how to evaluate, but don't have code to differentiate.  This is often because it's difficult/error-prone to write or because the interface by which we call it does not support derivatives.  (Commercial packages often fall in this category.)
Step11:  Automatically choosing a suitable $\epsilon$
Step12:  This algorithm is imperfect, leaving some scaling responsibility to the user.
It is the default in PETSc's "matrix-free" Newton-type solvers.
Tinkering in class (2016-09-06)
Step13:  In class 2016-09-14
We are looking for a function for which the linear fit (newton solver) is more reliable than a quadratic fit (cubic solver).
<ASSISTANT_TASK:>
Python Code:
%matplotlib notebook
from matplotlib import pyplot
import numpy
tests = []
@tests.append
def f0(x):
    return x*x - 2, 2*x
@tests.append
def f1(x):
    return numpy.cos(x) - x, -numpy.sin(x) - 1
@tests.append
def f2(x):
    return numpy.exp(-numpy.abs(x)) + numpy.sin(x), numpy.exp(-numpy.abs(x))*(-numpy.sign(x)) + numpy.cos(x)
@tests.append
def f3(x):
    return x*x - x + 0.25, 2*x - 1
@tests.append
def f4(x):
    return numpy.tan(x+2), numpy.cos(x+2)**(-2)
    return numpy.exp(-x*x), numpy.exp(-x*x)*(-2*x)
@tests.append
def f5(x):
    return (x <= 1)*1.0, 0*x
@tests.append
def f6(x):
    return x*numpy.sin(5/x), numpy.sin(5/x) - numpy.cos(5/x)*5/x
x = numpy.linspace(-2,2,100)
pyplot.plot(x, 0*x, color='k')
for f in tests:
    pyplot.plot(x, f(x)[0], label=f.__name__)
pyplot.legend(loc='upper right')
pyplot.style.use('ggplot')
pyplot.show()
Explanation: Jupyter notebooks
This is a Jupyter notebook using Python.  You can install Jupyter locally to edit and interact with this notebook.
Rootfinding
Rootfinding is the process of solving $$f(x) = 0$$ for $x$.  The standard assumption is that $f : R \to R$ is continuous.  We are interested in developing general-purpose algorithms---those that can use $f(x)$ as a black box without needing to look inside.  When we implement our rootfinding algorithm in software, the user will pass a function or program to compute $f(x)$.  Rootfinding methods for differentiable functions may also use the derivative $f'(x)$.
Some questions immediately arise:
* Existence. When does this equation have at least one solution?
* Uniqueness. When is the solution unique?
Let's consider some test functions, defined here along with their derivatives which we'll use later.
End of explanation
def hasroot(f, a, b):
    return f(a)[0]*f(b)[0] < 0
    
def bisect(f, a, b, verbose=False):
    mid = (a + b)/2.
    if b-a < 1e-5:
        return mid
    if verbose:
        print('bisect', mid)
    if hasroot(f, a, mid):
        return bisect(f, a, mid, verbose)
    else:
        return bisect(f, mid, b, verbose)
Explanation: Which of these functions have at least one root?
Which have more than one root?
Can we determine these properties merely by evaluating $f(x)$ for some values of $x$?
Bisection
Bisection is a rootfinding technique that starts with an interval $[a,b]$ containing a root and does not require derivatives.
End of explanation
bisect(tests[0], 0, 2)
numpy.sqrt(2) - bisect(tests[0], 0, 2)
Explanation: Notice that we need to define hasroot above.
Let's try running it:
End of explanation
bisect(tests[0], 0, 2, verbose=True)
Explanation: We get about 5 digits of accuracy.  Why?  How fast did we get there?
End of explanation
for f in tests:
    print(f.__name__, bisect(f, -2, 2.1))
Explanation: Can you find any problems with this implementation?  List them below:
No error checking
Can "converge" to nonsense
Recursion is bad in Python
Doesn't handle ill-defined points (function value not available)
Evaluates $f(x)$ more than necessary (multiple times at the same point)
Let's try running it on the rest of the test problem set:
End of explanation
def newton(f, x, verbose=False):
    for i in range(100):
        fx, dfx = f(x)
        if verbose:
            print(f.__name__, i, x, fx)
        if numpy.abs(fx) < 1e-12:
            return x, fx, i
        try:
            x -= fx / dfx
        except ZeroDivisionError:
            return x, numpy.NaN, i
for f in tests:
    print(f.__name__, newton(f, 1))
Explanation: What's going wrong here?  How can we improve the implementation and what are fundamental limitations of the algorithm?
Convergence rate
Let's quantitatively revisit the convergence rate.  A convergent rootfinding algorithm produces a sequence of approximations $x_i$ such that $$\lim_{i \to \infty} x_i \to x_$$ where $f(x_) = 0$.  For analysis, it is convenient to define the errors $e_i = x_i - x_$. We say that an algorithm is linearly convergent* if $$\lim_{i \to \infty} |e_{i+1}| / |e_i| = \rho < 1.$$  A smaller convergence factor $\rho$ represents faster convergence.
What is $\rho$ for bisection?
Remarks on bisection
Specifying an interval is often inconvenient
An interval in which the function changes sign guarantees convergence (robustness)
No derivative information is required
Roots of even degree are problematic
The solution error is directly available
The convergence rate is modest -- one iteration per bit of accuracy
Newton-Raphson Method
Much of numerical analysis reduces to Taylor series, the approximation
$$ f(x) = f(x_0) + f'(x_0) (x-x_0) + f''(x_0) (x - x_0)^2 / 2 + \dotsb $$
centered on some reference point $x_0$.
In numerical computation, it is exceedingly rare to look beyond the first-order approximation
$$ \tilde f_{x_0}(x) = f(x_0) + f'(x_0)(x - x_0) . $$
Since $\tilde f_{x_0}(x)$ is a linear function, we can explicitly compute the unique solution of $\tilde f_{x_0}(x) = 0$ as
$$ x = x_0 - f(x_0) / f'(x_0) . $$
This is Newton's Method (aka Newton-Raphson or Newton-Raphson-Simpson) for finding the roots of differentiable functions.
End of explanation
for f in tests:
    print(f.__name__, '{0:15.12f} {1:8.2e} {2:2d}'.format(*newton(f, -0.1)))
Explanation: Oops, how can we fix this?
This output is kinda hard to read, so let's make it cleaner.
End of explanation
def fquartic(x):
    return (x - 0.9)**4, 4*(x - 0.9)**3
newton(fquartic, 0)
Explanation: Did we solve all of these equations?
How can the iteration break down?
Does choosing a different initial guess lead to different solutions?
How is this convergence test different from the one we used for bisection?
Is the convergence rate similar for all test equations?
Convergence of Newton-type algorithms
We would like to know sharp conditions on when Newton-type algorithms converge, and if so, how fast.  This theory will build on that for a general Fixed Point Iteration $x_{i+1} = g(x_i)$ where $g$ is a continuously differentiable function.  Suppose that there exists a fixed point $r = g(r)$.  By the mean value theorem, we have that
$$ x_{i+1} - r = g(x_i) - g(r) = g'(c_i) (x_i - r) $$
for some $c$ with $|c - r| < |x_i - r|$.
In other words, $|e_{i+1}| = |g'(c_i)| |e_i|$, which converges to zero if $|g'(c_i)| < 1$.
If $|g'(r)| < 1$ then for any $\epsilon > 0$ there is a neighborhood of $r$ such that $|g'(c)| < |g'(r)| + \epsilon$ for all $c$ in that neighborhood.
Consequently, we have:
Theorem (Sauer 1.6): Linear Convergence of Fixed Point Iteration
If $g$ is continuously differentiable, $r = g(r)$, and $|g'(r)| < 1$ then the fixed point iteration $x_{i+1} = g(x_i)$ is locally linearly convergent with rate $|g'(r)|$.
Observations
A rootfinding problem $f(x) = 0$ can be converted to a fixed point problem $x = x - f(x) =: g(x)$ but there is no guarantee that $g'(r) = 1 - f'(r)$ will have magnitude less than 1.
Problem-specific algebraic manipulation can be used to make $|g'(r)|$ small.
$x = x - h(x)f(x)$ is also a valid formulation for any $h(x)$ bounded away from $0$.
Can we choose $h(x)$ such that $ 1 - h'(x)f(x) - h(x)f'(x) = 0$ when $f(x) = 0$?
In other words,
$$ x_{i+1} = x_i - ??? . $$
It turns out that Newton's method has locally quadratic convergence to simple roots, $\lim_{i \to \infty} |e_{i+1}|/|e_i^2| < \infty$.
"The number of correct digits doubles each iteration."
Now that we know how to make a good guess accurate, the effort lies in getting a good guess.
Culture: fast inverse square root
The following code appeared literally (including comments) in the Quake III Arena source code (late 1990s).
```C
float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
x2 = number * 0.5F;
y  = number;
i  = * ( long * ) &amp;y;                       // evil floating point bit level hacking
i  = 0x5f3759df - ( i &gt;&gt; 1 );               // what the fuck? 
y  = * ( float * ) &amp;i;
y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
return y;
}
```
We now have vector instructions for approximate inverse square root.
More at https://en.wikipedia.org/wiki/Fast_inverse_square_root
Conditioning
End of explanation
format((.2 - 1/3) + 2/15, '.20f')
# format((.2 - 1/3) + (1/3 - 0.2), '.20f')
# format((1 + 1e-12) - 1, '.20f')
eps = 1
while 1 + eps > 1:
    eps /= 2
eps_machine = eps  # We call this "machine epsilon"
numpy.log(1 + 1e-12) - numpy.log1p(1e-12)
(numpy.log(1 + 1e-12) - numpy.log1p(1e-12)) / numpy.log1p(1e-12)
x = numpy.array([1,1e5,1e10,1e15])
numpy.sin(numpy.pi*x)
numpy.sin(x)**2 + numpy.cos(x)**2 - 1
[numpy.tan((3.14159+eps)/2) for eps in [1e-6,1e-8]], 1/numpy.cos(3.14159)**2
Explanation: We only get three digits correct despite a very small residual (and it takes many iterations to get there).
Difficulty computing zeros of polynomials can also arise when all the roots are simple.  For example, the Wilkinson polynomial
$$ \prod_{i=1}^{20} (x - i) = \sum_{i=0}^{20} a_i x^i $$
has roots at each of the positive integers up to 20, but the roots are extremely sensitive to perturbations of the coefficients $a_i$, as shown in this figure from Trefethen and Bau (1999).
Numerical difficulties in which "correct" algorithms produce unreliable solutions almost always stem from lack of stability and/or ill conditioning.
Absolute condition number
Consider a function $f: X \to Y$ and define the absolute condition number
$$ \hat\kappa = \lim_{\delta \to 0} \max_{|\delta x| < \delta} \frac{|f(x + \delta x) - f(x)|}{|\delta x|} = \max_{\delta x} \frac{|\delta f|}{|\delta x|}. $$
If $f$ is differentiable, then $\hat\kappa = |f'(x)|$.
Floating point arithmetic
Floating point arithmetic $x \circledast y := \text{float}(x * y)$ is exact within a relative accuracy $\epsilon_{\text{machine}}$.  Formally,
$$ x \circledast y = (x * y) (1 + \epsilon) $$
for some $|\epsilon| \le \epsilon_{\text{machine}}$.
End of explanation
def diff(f, x, epsilon=1e-5):
    return (f(x + epsilon) - f(x)) / epsilon
diff(numpy.sin, 0.7, 1e-8) - numpy.cos(0.7)
x = .5
diff(numpy.tan, x) - 1/numpy.cos(x)**2
x = 3.14/2
[(eps, diff(numpy.tan, x, eps) - 1/numpy.cos(x)**2) for eps in [1e-14, 1e-12, 1e-10, 1e-8, 1e-6, 1e-4]]
x = 1e4
[(eps, diff(numpy.log, x, eps) - 1/x) for eps in [1e-14, 1e-12, 1e-10, 1e-8, 1e-6, 1e-4, 1e-2]]
Explanation: Relative condition number
Given the relative nature of floating point arithmetic, it is more useful to discuss relative condition number,
$$ \kappa = \max_{\delta x} \frac{|\delta f|/|f|}{|\delta x|/|x|}
= \max_{\delta x} \Big[ \frac{|\delta f|/|\delta x|}{|f| / |x|} \Big] $$
or, if $f$ is differentiable,
$$ \kappa = \max_{\delta x} |f'(x)| \frac{|x|}{|f|} . $$
How does a condition number get big?
Take-home message
The relative accuracy of the best-case algorithm will not be reliably better than $\epsilon_{\text{machine}}$ times the condition number.
$$ \max_{\delta x} \frac{|\delta f|}{|f|} \ge \kappa \cdot \epsilon_{\text{machine}} $$
Numerical differentiation
Suppose we want to apply Newton's method to a function that we know how to evaluate, but don't have code to differentiate.  This is often because it's difficult/error-prone to write or because the interface by which we call it does not support derivatives.  (Commercial packages often fall in this category.)
End of explanation
def diff_wp(f, x, eps=1e-8):
    Numerical derivative with Walker and Pernice (1998) choice of step
    h = eps * (1 + abs(x))
    return (f(x+h) - f(x)) / h
x = 1
[(eps, diff_wp(numpy.log, x, eps) - 1/x) for eps in [1e-14, 1e-12, 1e-10, 1e-8, 1e-6, 1e-4, 1e-2]]
x = 1e-4
[(eps, diff_wp(numpy.log, x, eps) - 1/x) for eps in [1e-14, 1e-12, 1e-10, 1e-8, 1e-6, 1e-4, 1e-2]]
Explanation: Automatically choosing a suitable $\epsilon$
End of explanation
numpy.log(-1)
x = numpy.sqrt(-1)
1/numpy.inf
x = numpy.linspace(0,3,100)
pyplot.plot(x, numpy.sqrt(x + 1e-1))
pyplot.show()
numpy.tan(1e100)
numpy.tan(1e100*(1 + 2*eps_machine))
1e100 - 1e100*(1 + 2*eps_machine)
Explanation: This algorithm is imperfect, leaving some scaling responsibility to the user.
It is the default in PETSc's "matrix-free" Newton-type solvers.
Tinkering in class (2016-09-06)
End of explanation
x = numpy.linspace(-5,5)
pyplot.figure()
f = numpy.sin(x) + 1.1*x
pyplot.plot(x, f)
Explanation: In class 2016-09-14
We are looking for a function for which the linear fit (newton solver) is more reliable than a quadratic fit (cubic solver).
End of explanation
<END_TASK>