Cours ENPC — Pratique du calcul scientifique
$$
% % %
%
$$
Goal of this chapter
Study numerical methods to find solutions to \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{0}}, where \mathbf{\boldsymbol{f}}\colon \mathbf R^n \to \mathbf R^n (number of unknows = number of equations).
Example applications
Calculate \sqrt{2} by solving x^2 = 2;
Solve two point boundary value problems, e.g.
u''(t) = f\bigl(t, u(t)\bigr), \qquad u(0) = u(1) = 0.
Solve nonlinear partial differential equations, e.g. \nabla \cdot \bigl(\alpha(u)\nabla u \bigr) = f.
Remark
A linear system of the form \mathsf A \mathbf{\boldsymbol{x}} = \mathbf{\boldsymbol{b}} can be recast in this form with \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}) = \mathsf A \mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{b}}. With this notation, Richardson’s method for linear systems rewrites \mathbf{\boldsymbol{x}}_{k+1} = \mathbf{\boldsymbol{x}}_k - \omega \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}_k). Questions:
Is such a method convergent for general nonlinear equations?
How fast is the convergence?
What is the optimal value of \omega.
Setting and idea
Suppose that f\colon [a, b] \to \mathbf R is continuous and f(a) \leqslant 0 \leqslant f(b).
{\color{green} \rightarrow} Then there is a root of f in [a, b]. Let x = \frac{a + b}{2}.
If f(a)f(x) \leqslant 0, then there is a root of f in [a, x];
If f(a)f(x) > 0, then there is a root of f in [x, b].
Convergence
It holds that |b_i - a_i| = \frac{|b - a|}{2^i} Both sequences (a_i) and (b_i) converge to a root x_* of f.
~
Comments
Works only in the one-dimensional setting;
Robust but slow…
Definition
We say that a sequence (\mathbf{\boldsymbol{x}}_k)_{k\in \mathbf N} converges to \mathbf{\boldsymbol{x}}_*
linearly if \lim_{k\to\infty} \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} \in (0, 1).
superlinearly if \lim_{k\to\infty} \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} = 0.
with order q if \mathbf{\boldsymbol{x}}_k \to \mathbf{\boldsymbol{x}}_* and \lim_{k\to\infty} \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}^q} \in (0, \infty). If q = 2, we say the sequence converges quadratically.
Examples
Most iterative methods for linear equations converge linearly.
The bisection method converges linearly.
Idea
Rewrite
\mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{0}} \qquad \Leftrightarrow \qquad \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}}
and use iterative scheme
\mathbf{\boldsymbol{x}}_{k+1} = \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_k)
Examples (more on these in the next section)
Chord method: for for \alpha \neq 0 \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} - \frac{\mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}})}{\alpha} {\color{green} \rightarrow} In the linear setting, this is Richardson’s iteration with \omega = \alpha^{-1}.
Newton–Raphson method: \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} - \mathsf J_f(\mathbf{\boldsymbol{x}})^{-1} \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}), where \mathsf J_f(\mathbf{\boldsymbol{x}}) is the Jacobian matrix of \mathbf{\boldsymbol{f}} at \mathbf{\boldsymbol{x}}.
Proposition: convergence under Lipschitz assumption
If \mathbf{\boldsymbol{F}} is Lipschitz continuous with L < 1, \forall (\mathbf{\boldsymbol{x}}, \mathbf{\boldsymbol{y}}) \in \mathbf R^n \times \mathbf R^n, \qquad {\lVert {\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{y}})} \rVert} \leqslant L {\lVert {\mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{y}}} \rVert} then
There exists a unique fixed point \mathbf{\boldsymbol{x}}_* of \mathbf{\boldsymbol{F}}.
The sequence (\mathbf{\boldsymbol{x}}_k)_{k\in \mathbf N} converges exponentially to \mathbf{\boldsymbol{x}}_*: \forall k \in \mathbf N, \qquad {\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert} \leqslant L^k {\lVert {\mathbf{\boldsymbol{x}}_0 - \mathbf{\boldsymbol{x}}_*} \rVert}.
Sketch of proof
Show that (\mathbf{\boldsymbol{x}}_k) is Cauchy, thus convergent to some \mathbf{\boldsymbol{x}}_* \in \mathbf R^n.
Show that \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*) = \mathbf{\boldsymbol{x}}_*.
Deduce exponential convergence from \mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_* = \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_k) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*).
Proposition: local convergence under local Lipschitz condition
Assume that \mathbf{\boldsymbol{x}}_* is a fixed point of \mathbf{\boldsymbol{F}} and \forall \mathbf{\boldsymbol{x}} \in B_{\delta}(\mathbf{\boldsymbol{x}}_*), \qquad {\lVert {\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*)} \rVert} \leqslant L {\lVert {\mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{x}}_*} \rVert} \qquad(1) for L < 1. If \mathbf{\boldsymbol{x}}_0 \in B_{\delta}(\mathbf{\boldsymbol{x}}_*), then \mathbf{\boldsymbol{x}}_k \in B_{\delta}(\mathbf{\boldsymbol{x}}_*) for all k \in \mathbf N and {\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert} \leqslant L^k {\lVert {\mathbf{\boldsymbol{x}}_0 - \mathbf{\boldsymbol{x}}_*} \rVert}.
~
Proposition: local convergence criteria based on the Jacobian
Assume that \mathbf{\boldsymbol{F}} is differentiable at a fixed point \mathbf{\boldsymbol{x}}_*.
There is \delta > 0 such that Equation 1 is satisfied iff {\lVert {\mathsf J_F(\mathbf{\boldsymbol{x}}_*)} \rVert} < 1.
There is \delta > 0 such that Equation 1 is satisfied in some induced norm iff \rho\bigl(\mathsf J_F(\mathbf{\boldsymbol{x}}_*)\bigr) < 1.
Proposition: superlinear convergence
Assume that \mathbf{\boldsymbol{x}}_* is a fixed point of \mathbf{\boldsymbol{F}} and that \mathsf J_F(\mathbf{\boldsymbol{x}}_*) = 0. If \mathbf{\boldsymbol{x}}_k \to \mathbf{\boldsymbol{x}}_* as k \to \infty, then \lim_{k \to \infty} \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} = 0.
Remark
Under the assumption, there is \delta > 0 such that (\mathbf{\boldsymbol{x}}_k)_{k\geqslant 0} is a sequence converging to \mathbf{\boldsymbol{x}}_* for all \mathbf{\boldsymbol{x}}_0 \in B_{\delta}(\mathbf{\boldsymbol{x}}_*).
Proof
It holds that \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} = \frac{{\lVert {\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_{k}) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*)} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} = \frac{{\lVert {\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_{k}) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*) - \mathsf J_F(\mathbf{\boldsymbol{x}}_*) (\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*)} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}}. Since \mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_* \to \mathbf{\boldsymbol{0}} as k \to \infty, the right-hand side converges to 0 by definition of the Jacobian matrix.
Description of the method
The chord method is the iteration \mathbf{\boldsymbol{x}}_{k+1} = \mathbf{\boldsymbol{x}}_k - \mathsf A^{-1} \mathbf{\boldsymbol{f}}(x_k), with \mathsf A \in \mathbf R^{n \times n} an invertible matrix parameter.
~
Remarks
This is a fixed-point method with \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} - \mathsf A^{-1} \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}).
Clearly \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} if and only if \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{0}}.
When n = 1, the iteration reads x_{k+1} = x_k - \frac{f(x_k)}{\alpha}, {\color{green} \rightarrow} In the linear setting, this is Richardson’s iteration.
Convergence study when n=1
In this case F(x) = x - \frac{f(x)}{\alpha}, \qquad F'(x) = 1 - \frac{f'(x)}{\alpha} We deduce:
The iteration converges locally close to the fixed point iff \left| 1 - \frac{f'(x_*)}{\alpha} \right| < 1. {\color{green} \rightarrow}\alpha must be of the same sign as f'(x_*) and sufficiently large: \left|\alpha\right|>\frac{\left|f'(x_*)\right|}{2}
The convergence is superlinear if \alpha = f'(x_*)
F(x) = x - \frac{1}{\alpha}\tanh(x-1), \qquad x_* = 1,
F'(x_*) = 1 - \frac{1}{\alpha\,\cosh^2(x-1)}=1-\frac{1}{\alpha}
F(x) = x - \frac{1}{\alpha}\tanh(x-1), \qquad x_* = 1,
F'(x_*) = 1 - \frac{1}{\alpha\,\cosh^2(x-1)}=1-\frac{1}{\alpha}
F(x) = x - \frac{1}{\alpha}\tanh(x-1), \qquad x_* = 1,
F'(x_*) = 1 - \frac{1}{\alpha\,\cosh^2(x-1)}=1-\frac{1}{\alpha}
F(x) = x - \frac{1}{\alpha}\tanh(x-1), \qquad x_* = 1,
F'(x_*) = 1 - \frac{1}{\alpha\,\cosh^2(x-1)}=1-\frac{1}{\alpha}
Geometric interpretation
At each iteration:
Approximate f(x) by the straight line \widehat f_k(x) = f(x_k) + \alpha (x - x_k)
Let x_{k+1} be the root of \widehat f_k.
Toy problem
Let us apply the chord method to f(x) = x^2 - 2. In this case F(x) = x - \frac{x^2 - 2}{\alpha}, \qquad F'(x) = 1 - \frac{2x}{\alpha}. Here x_* = \sqrt{2} and f'(x_*) = 2\sqrt{2}.
Motivation
The chord method is optimal with \mathsf A = \mathsf J_f(\mathbf{\boldsymbol{x}}_*);
However, \mathbf{\boldsymbol{x}}_* is unknown a priori…
Newton-Raphson method: use \mathsf A = \mathsf A(k) = \mathsf J_f(\mathbf{\boldsymbol{x}}_k). which leads to the iteration \mathbf{\boldsymbol{x}}_{k+1} = \mathbf{\boldsymbol{x}}_k - \mathsf J_f(\mathbf{\boldsymbol{x}}_k)^{-1} \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}_k)
One-dimensional version
In dimension one, the Newton-Raphson iteration reads x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}
Remarks
The method requires to solve a linear system at each step.
The iteration is well defined only if \mathbf{\boldsymbol{J}}_f(\mathbf{\boldsymbol{x}}_k) is nonsingular.
The method is a fixed point iteration with \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} - \mathsf J_f(\mathbf{\boldsymbol{x}})^{-1} \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}).
Dimension 1: If f'(x_*) \neq 0, then F'(x_*) = 0. {\color{green} \rightarrow} local superlinear convergence.
Dimension n: If \mathsf J_f(\mathbf{\boldsymbol{x}}_*) is nonsingular, then \mathsf J_F(\mathbf{\boldsymbol{x}}_*) = \mathsf 0. {\color{green} \rightarrow} local superlinear convergence.
Theorem: Quadratic convergence
Assume that f \in C^2(\mathbf R), that f'(x_*) \neq 0, and that the Newton-Raphson iteration converges to x_*. Then \lim_{k \to \infty} \frac{{\lvert {x_{k+1} - x_*} \rvert}}{{\lvert {x_{k} - x_*} \rvert}^2} = \left| \frac{f''(x_*)}{2 f'(x_*)} \right|.
Proof
We note that x_{k+1} - x_* = x_{k} - \frac{f(x_k)}{f'(x_k)} - x_* = \frac{1}{f'(x_k)} \underbrace{\Bigl(f'(x_k)(x_{k} - x_*) - f(x_k)\Bigr)}_{= \frac{1}{2} f''(\xi_k) (x_{k} - x_*)^2}. for some \xi_k between x_k and x_*. We conclude that x_{k+1} - x_* = \frac{f''(\xi_k) (x_k - x_*)^2}{2f'(x_k)}.
Toy problem
Let us apply the chord method to f(x) = (x^2 - 2)(x - 3)^4.
Single roots \sqrt{2} and -\sqrt{2}.
Quadruple root 3.
{\color{green} \rightarrow} superlinear convergence does not apply.
~
~
~
Geometric interpretation
At each iteration:
Approximate f(x) by its tangent at x_k \widehat f_k(x) = f(x_k) + f'(x_k) (x - x_k)
Let x_{k+1} be the root of \widehat f_k.
