THREE FUNDAMENTAL LEVELS OF NUMERICAL MATHEMATICS

Numerical mathematics ≈ practical (physical, machinane aided) evaluation of theoretical formulas

1) Pure mathematical level Example. Riemann conjecture ζ(z) = ∑_{n=1}^{∞} \frac{1}{n^{z}}

ζ(z) = 0, z ≠ −2, −4, … ⟹ ∃ t ∈ IR z = \frac{1}{2} + it

Compromise-free axionmatic poof is required. Alan TURING: looking for counter-examples with (ancient) computer.

2) Theoretical analysis of computation algorithms.

Theoretical accuracy, reliability of guesses for solutions. (”Usual” university texts)

Investigation of finite computations with ”Infinitely long arithmetics”. Example. Convergence rate estimates for Newton iteration.

We are going farther: investigation of possible chaos in calculations.

3) ”Engineering” level.

Use of machine architectures and algorithms without theoretical criticism.

Not for disprising even from the view points of 1-2):

Learning- and genetic algorithms — not yet controlled, but extremally successful. Topics for 2). Quantum computers

Example for the cooperation of 1-2-3).

Computations in Quantum Chemistry (Gaussian package, Nobel Prize 2000).

1) Schrödinger-equation,

2) Approximating solutions techniques (⟶ conjugate gradient methods).

3) Imitation with computer of 2) (far from financial success in the moment)

BASIC CONCEPTS IN NUMERICAL MATHEMATICS

Task: F (x, d) = 0 where d := [data] Direct problem: given F, d, x =? Inverse problem: given F, x, d ∈ {?}; Identification problem: given x, d, F =?. Well-posedness: existing unique solution with continuous dependence on the data Relative condition number (for F (x + δx, d + δd) = 0): K(d) := sup_{δd} \frac{∥δx∥/∥x∥}{∥δd∥/∥d∥}. Condition number (when d = 0 or x = 0): K_{abs}(d) := sup_{δd} \frac{∥δx∥}{∥δd∥}. Stability for approximating methods with F_{n}(x_{n}, d_{n}) = 0:

F_{n}(x, d) → 0 (n → ∞) whenever F (x, d) = 0.

A proiori analysis: Investigation of direct problems. A posteriori analysis: Investigation of indirect problems.

NEWTON ITERATION

Recall the classical 1D case. Mesopotamia (≈ 1500 BC). \sqrt{2}: 1+\frac{30}{60} =1.5 1+\frac{25}{60} ≈1.41666 1+\frac{24}{60}+ \frac{51}{3600} ≈1.1.414667 2, \frac{1}{2}[2 + \frac{2}{2}], \frac{1}{2}[1.S + \frac{2}{1.5}], \frac{1}{2}[1.41˙6 + \frac{2}{1.41˙6}]

x^{2} = 2 x = \frac{2}{x} \begin{matrix}x = \frac{2}{x} \\ x = x \end{matrix} \}·1/2_{\frac{·1/2}{+}} x = \frac{1}{2}(x + \frac{2}{x})

\begin{matrix}\sqrt{a} − ra − is \\ x^{2} = a \end{matrix} \begin{matrix}x = \frac{a}{x} \\ x = x \end{matrix} \} x = \frac{1}{2}[x + \frac{a}{x}]

Questions

1) Why does it converge

2) How to proceed with other numbers, e.g. ^{3} \sqrt{a} =?

KEPLER \begin{matrix}x^{3} = a \\ x = x \end{matrix} \}· \frac{2}{3}^{\frac{1}{3}} x = \frac{2}{3}x + \frac{1}{3} \frac{a}{x^{2}}

Conjecture ^{N} \sqrt{a}-hz x_{0} := a > 1 x_{n+1} = (1 − \frac{1}{N} )x_{n} + \frac{1}{N} \frac{a}{x_{n}^{N−1}}

1) NEWTON’s answer: in general for the solution of f(x) = 0 [in particular x^{2} − 2 = 0 ITERATION WITH TANGENT LINES Slope = f′(x_{n}) FIGURE (x_{n+1} − x_{n}) · f′(x_{n}) = f(x_{n})

x_{n+1} = x_{n} − \frac{f(x_{n})}{f′(x_{n})}

Example. ^{N} \sqrt{a} x^{N} − a = 0 f(x) := x^{N} − a f′(x) = Nx^{N−1}

x_{n+1} = x_{n} − \frac{x_{n}^{N} − a}{Nx_{n}^{N−1}} = (1 − \frac{1}{N} )x_{n} + \frac{1}{N} \frac{a}{x_{n}^{N−1}}

Reasons for convergence? Not always feasable: f(x) := x^{3} − x − 1/\sqrt{3} ”finding 0” may be CHAOTIC [details later] TAYLOR FORMULA + LAGRANGE REMAINDER TERM

\begin{aligned}f(x + h) & = f′(x) + f′(x)h + \frac{1}{2}f′′(x)h^{2} + … + \frac{1}{(d − 1)!}f^{(d−1)}(x)h^{d−1}+ \\ & + \frac{1}{d!}f^{(d)}(t_{x,h})h^{d} ∃ t_{x,h} ∈ [x, x + h] , \end{aligned}

if f : IR→IR is d times continuously differentiable. d = 2: f(x + h) = f(x) + f′(x)h + \frac{1}{2}f′′(t_{x,h})h^{2} NEWTON’s HEURISTICS x_{n} ≈ x_{∗} f(x_{∗}) = 0 f(x_{n} + h) ≈ f(x_{n}) + f′(x_{n})h 0 = f(x_{∗}) ≈ f(x_{n}) + f′(x_{n}) · (x_{∗} − x_{n}) x_{∗} ≈ x_{n} − f(x_{n})/f′(x_{n}) ESTIMATE

(1) 0 = f(x_{∗}) = f(x_{n}) + f′(x_{n})(x_{∗} − x_{n}) + \frac{1}{2}f′′(t_{x_{n},x_{n}−x_{∗}})(x_{∗} − x_{n})^{2}

(2) 0 = [f replaced with its linear approximation ](x_{n+1}) = f(x_{n}) + f′(x_{n})(x_{n+1} − x_{n})

(2) - (1) ⇒ 0 = f′(x_{n})(x_{n+1} − x_{∗}) − \frac{1}{2}f′′(t_{x_{n},x_{∗}−x_{n}})(x_{n} − x_{∗})^{2}

|x_{n+1} − x_{∗}| = \frac{1}{2} \underbrace{\frac{|f′′(t_{x_{n},x_{n}−x_{∗}})|}{|f′(x_{n})|}}_{≤ M , ha x_{n} ∈ (x_{∗} − ε, x_{∗} + ε).} |x_{n} − x_{∗}|^{2}

h_{n} := |x_{n} − x_{∗}| = [ distance of x_{n} from the solution ] h_{n+1} ≤ Mh_{n}^{2} , but M can be very large Iterating this formally yields

h_{1} ≤ Mh_{0}^{2}
h_{2} ≤ Mh_{1}^{2} ≤ M(Mh_{0}^{2}) = M^{3}h_{0}^{4}
h_{3} ≤ Mh2_{2} ≤ M(M^{3}h_{0}^{4}) = M^{7}h_{0}^{8}

⋮

h_{n+1} ≤ M^{2^{n}−1}h_{0}^{2^{n}} = \frac{1}{M} [Mh_{0}]^{2^{n}}

If Mh_{0} < 1 (x_{0} lies ”very near” to x_{∗}) then x_{n} → x_{∗}, moreover |x_{n} −x_{∗}| ≤ const.(1−δ)^{2^{n}}

NEWTON ITERATION IN SEVERAL VARIABLES

F : IR^{N} → IR^{N}, x_{i} : [ ξ_{1} ; ⋮ ; ξ_{N} ] ↦ ξ_{i} coordinates (coordinate functions), F = [ f_{1} ; ⋮ ; f_{N} ] F ′_{H}(A) := lim_{τ→0} \frac{1}{τ} [F (A + τH) − F (A)] = \frac{∂F (A + τH)}{∂τ}|_{τ=0} directional derivative Recall. F continuously differentiable ⇒ H ↦ F ′_{H}(A) is a linear mapping IR^{N} →IR^{N} F d times continuously differentiable ⇒ (H_{1}, …, H_{d}) ↦ F ′_{H_{1}} _{···}^{...}′_{H_{d}}(A) is a symmetric d-linear [IR^{N}]^{d} →IR^{N} mapping

Notation: F^{(d)}(A)H_{1} ⋯ H_{d} := F ′_{H_{1}} _{···}^{...}′_{H_{d}}(A), F^{(d)}(A)H^{d} := F^{(d)}(A) \underbrace{H ⋯ H}_{d times}

Tensorial form: F^{(d)} ∼ [\frac{∂^{d}F}{∂x_{k_{1}} ⋯ ∂x_{k_{d}}} ]_{k_{1},...,k_{d}=1}^{N}.

Taylor formulas with remainder terms

F_{[A,A+H]}^{(d)} := d! ∫_{t_{1}=0}^{1} ∫_{t_{2}=0}^{1} ⋯ ∫_{t_{d}=0}^{1} F^{(d)}(A + t_{d}H) dt_{d} dt_{d−1} ⋯ dt_{1};

F (A + H) = ∑_{k=0}^{d−1} \frac{1}{k!}F^{(k)}(A)H^{k} + \frac{1}{d!}F_{[A,A+H]}^{(d)}H^{d} =

= ∑_{k=0}^{d−1} \frac{1}{k!}F^{(k)}(A)H^{k} + \frac{1}{d!} ∫_{τ=0}^{1} \underbrace{w_{d}(τ)}_{d(1−τ)^{d−1}} F^{(d)}(A + τH)H^{d} dτ, ∫_{0}^{1}w_{d} = 1;

= [F^{(d)}(A + ϑ_{1}H), …, F^{(d)}(A + ϑ_{N}H) convex lin. combination] =

= ∑_{k=0}^{d−1} \frac{1}{k!}F^{(k)}(A)H^{k} + \frac{1}{d!}⟨ϕ|F^{(d)}(A + ϑ_{A,H}H)H^{d}⟩U, 1 = ⟨ϕ|U⟩ = ∥ϕ∥ = ∥U∥.

Iteration step

X_{n+1} := X_{n} − [F′(X_{n})]^{−1}F (X_{n})

That is X_{n+1} is the solution of the equation [Taylor polynomial of 1st order of F around X_{n}] = 0:

(∗)  F (X_{n}) + F′(X_{n})(X_{n+1} − X_{n}) = 0 .

This is well6-defined and unique if the linear map F′(X_{n}) : IR^{N} → IR^{N} is invertible.

Convergence estimate: Similarly as with 1 variable, but using Taylor formula with remainder in integration:

0 = F (X_{n}) + F′(X_{n})(X_{n+1} − X_{n}),
0 = F (X_{∗}) =

= F (X_{n}) + F′(X_{n})(X_{∗}−X_{n}) + \frac{1}{2} ∫_{τ=0}^{1} w_{2}(τ)F′′(X_{n} + τ(X_{∗}−X_{n}))(X_{∗}−X_{n})^{2} dτ.

Taking the differnce of these equations we get

0 = F′(X_{n})(X_{n+1}−X_{∗}) − \frac{1}{2} ∫_{τ=0}^{1} w_{2}(τ)F′′(X_{n} + τ(X_{∗}−X_{n}))(X_{∗}−X_{n})^{2} dτ,

X_{n+1} − X_{∗} = \frac{1}{2} ∫_{τ=0}^{1} w_{2}(τ)[F′(X_{n})]^{−1}F′′(X_{n} + τ(X_{∗}−X_{n}))(X_{∗}−X_{n})^{2} dτ,

∥X_{n+1} − X_{∗}∥ ≤ \frac{1}{2} max_{τ∈[0,1]} ∥[F′(X_{n})]^{−1}F′′(X_{n} + τ(X_{∗}−X_{n}))(X_{∗}−X_{n})^{2}∥ ≤

≤ \frac{1}{2}[ max_{{X,Y ∈[Xn,X∗] \atop ∥H∥=1}} ∥[F′(X)]^{−1}F′′(Y)H^{2}∥]∥X_{∗} − X_{n}∥^{2}.

Theorem. Assume F : IR^{N} → IR^{N} is a 2 times continuously differntiable mapping, such that F (X_{∗}) = 0 at the point X_{∗} ∈IR^{N}. Suppose furthermore that K is a bounded sphericak neighvorhood of X_{∗} with the property that det(F′(X)) ≠ 0 (X ∈K). Then the constant

M := \frac{1}{2} max_{{X,Y ∈K \atop ∥H∥=1}} ∥[F′(X)]^{−1}F′′(Y)H^{2}∥

is a well-defines finite number. If X_{0} ∈ K lies so near to X_{∗} that we have

λ := M∥X_{0} − X_{∗}∥ < 1 ,

then all the points X_{1}, X_{2}, … are (uniquely) well-defined by the steps (∗) of the Newton iteration starting from X_{0} and they remain in K (i.e. X_{1}, X_{2}, … ∈ K) and ∥X_{n} − X_{∗}∥ = O(λ^{2^{n}}). Namely

∥X_{n} − X_{∗}∥ ≤ λ^{2^{n}−1}∥X_{0} − X_{∗}∥ (n = 0, 1, 2, …).

Proof. Verification that the constant M is well-defined: since F is 2 times continuously differentiable, the function (X, Y, H) ↦ ∥[F′(X)]^{−1}F′′(Y)H^{2}∥ is continuous and hence it assumes it maximum felveszi in the closed bounded figure K × K × {H : ∥H∥ = 1} ⊂ [IR^{N}]^{3}. Let K := {X : ∥X∥ ≤ ε}, X_{0} ∈ K és λ = M∥X_{0} − X_{∗}∥ < 1. Hnceforth we can use analogous arguments as in the case of 1 variable. We have ∥X_{0} − X_{∗}∥ = λ^{2^{0}−1}∥X_{0} − X_{∗}∥ ≤ ε. Induction step in accordance with the convegence estimate:

X_{n} ∈K + ∥X_{n}−X_{∗}∥ ≤ λ^{2^{n}−1}∥X_{0} − X_{∗}∥ ⟹
∥X_{n+1} − X_{∗}∥ ≤ M∥X_{n} − X_{∗}∥^{2} ≤ M[λ^{2^{n}−1}∥X_{0} − X_{∗}∥]^{2} =^{M∥X_{0}−X_{∗}∥=λ}

= λ^{2(2^{n}−1)+1}∥X_{0}−X_{∗}∥ = λ^{2^{n+1}−1}∥X_{0}−X_{∗}∥ ≤ ε (⇒ X_{n+1} ∈K).

Procedure with matrix calculus

∇f_{i}(x) = [\frac{∂f_{i}}{∂x_{1}} \frac{∂f_{i}}{∂x_{2}} … \frac{∂f_{i}}{∂x_{N}} ] gradient (row vector)

∇^{2}f_{i}(x) := [\frac{∂^{2}f_{i}}{∂x_{k}∂x_{ℓ}} ]_{k,ℓ=1}^{N} Hessian (matrix)

f_{i}(x + h) = f_{i}(x) + [∇f_{i}(x)]h⟩ + \frac{1}{2}h^{∗}[∇^{2}f_{i}(t_{x,h})]h

Suppose the point X_{n} ∈ IR^{N} is a given [ n-th approximation of X_{∗} where F (X_{∗}) = 0 ]

f_{i}(X) ≈ f_{i}^{[n]}(X) := [ approximationin1 − storderoff_{i} around X_{n}] =
= f_{i}(X_{n}) + [∇f_{i}(X_{n})](X − X_{n}) =
= f_{i}(x_{n}) + ⟨[∇f_{i}(x_{n})]^{∗}| X − X_{n}⟩

F (X) = 0 ≈ f_{i}^{[n]}(x_{1}, …, x_{n}) = 0 (i = 1, 2, …, N) system of equations

Example. N = 2, X = [ \begin{matrix}x \\ y \end{matrix} ], F = [ \begin{matrix}f \\ g \end{matrix} ], X = [ \begin{matrix}x_{n} \\ y_{n} \end{matrix} ]

\begin{matrix}F (X_{n}) + \frac{∂f(X_{n})}{∂x}(X_{n+1} − X_{n}) + \frac{∂f(X_{n})}{∂y}(y_{n}, −y_{n}) & = 0 \\ g(X_{n}) + \frac{∂g(X_{n})}{∂x}(x_{n+1} − x_{n}) + \frac{∂g(X_{n})}{∂y}(y_{n+1} − y_{n}) & = 0 \end{matrix} \}

f(X_{n}) + \underbrace{[ \begin{matrix}∇f \\ ∇g \end{matrix} ]}_{2 × 2 matrix the Jacobian of F} (X_{n+1} − X_{n}) = 0

In general: (in arbitrary dimensions N ≥ 1)

F′(X) = [ Jacobian of F at the point X] := \underbrace{[ ∇f_{1}(X) ; ⋮ ; ∇f_{N}(X) ]}_{N × N-es mátrix} .

We have to find the solution X_{n+1} of

F (X_{n}) + F′(X_{n})(X_{n+1} − X_{n}) = 0.

Hence

X_{n+1} := X_{n} − \underbrace{F′(X_{n})}_{Jacobi}^{−1}F (X_{n}) = X_{n} − [ ∇f_{1}(X_{n}) ; ⋮ ; ∇f_{N}(X_{n}) ]^{−1}F (X_{n})

Estimate. Since

f_{i}(X_{n} + H) = f_{i}(X_{n}) + [∇f_{i}(X_{n})]H + \frac{1}{2}H^{∗}[∇^{2}f_{i}(T_{i})]H ∃ T_{i} ∈ [X_{n}, X_{n} + H],

by setting H := X_{∗} − X_{n} we have

\begin{matrix}0 = f_{i}(X_{∗}) = f_{i}(X_{n}) + [∇f_{i}(X_{n})](X_{∗} − X_{n}) + \frac{1}{2}(X_{∗} − X_{N})^{∗}[∇^{2}f_{i}(T_{i})](X_{∗} − X_{N}), \\ 0 = F (X_{∗}) = F (X_{n}) + F′(X_{n})(X_{∗} − X_{n}) + \frac{1}{2} [ (X_{∗} − X_{N})^{∗}[∇^{2}f_{1}(T_{i})](X_{∗} − X_{N}) ; ⋮ ; (X_{∗} − X_{N})^{∗}[∇^{2}f_{1}(T_{i})](X_{∗} − X_{N}) ] \end{matrix}

for suitable T_{1}, …, T_{N} ∈ [X_{n}, X_{∗}]. On the other hand, by the definition of X_{n+1} we have

0 = F (X_{n}) + F′(X_{n})(X_{n+1} − X_{n}).

Subtracting the previous equations (vectorially) from this, we get

X_{n+1} − X_{∗} = \frac{1}{2}F′(X_{n})^{−1} [ (X_{n} − X_{∗})^{∗}∇^{2}f_{1}(T_{1})(X_{n} − X_{∗}) ; ⋮ ; (X_{n} − X_{∗})^{∗}∇^{2}f_{N}(T_{N})(X_{n} − X_{∗}) ] .

Theorem. If F is C^{2}-smooth, F (X_{∗}) = 0 and the Jacobian matrix F′(X_{∗}) is invertible then there is a (small) convex neighborhood U of the point X_{∗} and a (large) constant M such that starting from X_{0} ∈ U each step X_{n} is well-defined, belongs to U and

∥X_{n+1} − X_{∗}∥ ≤ M∥X_{n} − X_{∗}∥^{2}.

In particular

∥X_{n} − X_{∗}∥ ≤ μ^{2^{n}}/M → 0 if μ := M∥X_{0} − X_{∗}∥ < 1.

Example: 2D GPS-problem. Let S_{1}, S_{2}, S_{3} ∈ IR^{2} be three given distinct points. Determine the coordinates of a point P ∈ IR^{2} which satisfies the equiations

δ_{k} := d(P, S_{k}) − d(P, S_{k+1})) = ∥P − S_{k}∥ − ∥P − S_{k}∥ (k = 1, 2)

concerning distance differences.

Newton iteration for solving the equation F (P) = 0 where

F (X) := [ \begin{matrix}d(X, S_{1}) − d(X, S_{2}) − δ_{1} \\ d(X, S_{2}) − d(X, S_{3}) − δ_{2} \end{matrix} ] (X ∈ IR^{2}).

In terms of the canonical coordinates X ≡ (x, y), S_{k} ≡ (p_{k}, q_{k}) we can write

F (X) = [ \begin{matrix}\sqrt{(x − p_{1})^{2} + (y − q_{1})^{2}} − \sqrt{(x − p_{2})^{2} + (y − q_{2})^{2}} − δ_{1} \\ \sqrt{(x − p_{2})^{2} + (y − q_{2})^{2}} − \sqrt{(x − p_{3})^{2} + (y − q_{3})^{2}} − δ_{2} \end{matrix} ] .

Notice that the firsr order Taylor approximation of F around the point X_{n} ≡ (x_{n}, y_{n}) of the iteration is

\begin{aligned}F (X_{n} + V) & ≈^{(1)} F (X_{n}) + F′(X_{n})V = F (X) + \frac{d}{dτ} |_{τ=0}F (X_{n} + τV) = \\ & = F (x_{n}, y_{n}) + \frac{d}{dτ} |_{τ=0}F (x_{n} + τv_{n}, y_{n} + τw_{n}). \end{aligned}

In general,

\frac{d}{dτ} |_{τ=0}F (x + τv, y + τw) = [ \begin{matrix}\frac{(x−p_{1})v+(y−q_{1})w}{\sqrt{(x−p_{1})^{2}+(y−q_{1})^{2}}} − \frac{(x−p_{2})v+(y−q_{2})w}{\sqrt{(x−p_{2})^{2}+(y−q_{2})^{2}}} \\ \frac{(x−p_{2})v+(y−q_{2})w}{\sqrt{(x−p_{2})^{2}+(y−q_{2})^{2}}} − \frac{(x−p_{3})v+(y−q_{3})w}{\sqrt{(x−p_{3})^{2}+(y−q_{3})^{2}}} \end{matrix} ] .

If X_{n} is known then

X_{n+1} = X_{n} + V, where V is the solution of the linear equation F (X_{n}) + F′(X_{n})V = 0.

With given data: Let S_{1} ≡ (0, 10), S_{2} ≡ (8, 6), S_{3} ≡ (10, 0); assume P lies near to the point X_{0} ≡ (4.8, 6.4) (whose respective distances from S_{1}, S_{2}, S_{3} are 6, 3.225, 8.246), but δ_{1} = d(P, S_{1}) − d(P, S_{2}) = 3, δ_{2} = d(P, S_{2}) − d(P, S_{3}) = −5. Then

X_{1} ≡ (4.8 + v, 6.4 + w), where F (4.8, 6.4) + \frac{d}{dτ} |_{τ=0}F (4.8 + τu, 6 − 4 + τw) = 0,

Here we have

F (4.8, 6.4) = [ \begin{matrix}\sqrt{4.8^{2} + (6.4 − 10)^{2}} − \sqrt{(4.8 − 8)^{2} + (6.4 − 6)^{2}} − 3 \\ \sqrt{(4.8 − 8)^{2} + (6.4 − 6)^{2}} − \sqrt{(4.8 − 10)^{2} + 6.4^{2}} + 5 \end{matrix} ] = [ \begin{matrix}−0.225 \\ −0.021 \end{matrix} ],

\frac{d}{dτ} |_{τ=0}F (4.8 + τu, 6.4 + τw) = [ \begin{matrix}\frac{4.8v+(6.4−10)w}{6} − \frac{(4.8−8)v+(6.4−6)w}{3.225} \\ \frac{(4.8−8)v+(6.4−6)w}{3.225} − \frac{(4.8−10)v+6.4w}{8.246} \end{matrix} ].

Hence v = 0.092, w = −0.083, X_{1} ≡ (4.892, 6.317). Then d(X_{1}, S_{1}) = 6.123, d(X_{1}, S_{2}) = 3.124, d(X_{1}, S_{3}) = 8.124, The distance differences became improved: \overline{δ}_{1} = 2.999, \overline{δ}_{2} = −5.124.

Inner convergence estimate for Newton iteration D ⊂ IR^{N} compact region, f : D → IR^{N} C^{2}-smooth, f′(x) invertible (x ∈ D)

M := max \{\frac{1}{2}∥f′(x_{1})^{−1}f′′(x_{2})∥ : x_{1}, x_{2} ∈ D\}.

x_{0} ∈ D, x_{k−1} := x_{k} − f′(x_{k})f(x_{k}) well-defined (k = 1, …, n). Question. If x_{n} and x_{n+1} are ”close”, how far can be a root of f from x_{n+1}? Assumption: f(x_{∗}) = 0, [x_{n}, x^{∗}] ⊂ D

0 = f(x_{∗}) = f(x_{n}) + f′(x_{n})(x_{n} − x_{∗}) + \underbrace{\frac{1}{2}f′′_{[x_{n},x_{∗}]}}_{∫_{t_{1}=0}^{1} ∫_{t_{2}=0}^{1} f′′(x_{n}+t_{2}(x_{∗}−x_{n})) dt_{2} dt_{1}} (x_{∗} − x_{n})^{2}

y_{∗} := x_{n} − x_{∗}

y_{∗} = −f′(x_{n})^{−1}f(x_{n}) − \frac{1}{2}f′(x_{n})^{−1}f′′_{[x_{n},x_{∗}]}y_{∗}^{2}

y_{∗} = T (y_{∗}), T (y) := \underbrace{−f′(x_{n})^{−1}f(x_{n})}_{x_{n+1} − x_{n}} −\frac{1}{2}f′(x_{n})^{−1}f′′_{[x_{n},x_{∗}]}y^{2}

Recall. According to Brower’s Fixed Point Theorem, if a continuous mapping T maps a closed ball (or a topologically equivalent figure to a ball) into itself then T admits a fixed point (T (x) = x). Hence if for some δ > 0 we have T : \overline{B(δ)} = \{y ∈ IR^{N} : ∥y∥ ≤ δ\} → \overline{B(δ)}, then ∃ y_{∗} ∈ \overline{B(δ)} y_{∗} = T (y_{∗}).

Thus if there exists δ > 0, such that x_{n} + \overline{B(δ)} ⊂ D and ∥T (y)∥ ≤ δ whenever ∥y∥ ≤ δ, then, for any value δ with the above property, we can find a point x_{∗} ∈ D, such that f(x_{∗}) = 0 and ∥x_{n} − x_{∗}∥ ≤ δ. How large can be the value of such a δ?

∥T (y)∥ ≤ δ (∥y∥ ≤ δ), ⇐ ∥x_{n+1} − x_{n}∥ + Mδ^{2} ≤ δ,

Mδ^{2} − δ + ∥x_{n+1} − x_{n}∥ ≤ 0,

δ ∈ [\frac{1−\sqrt{1−4M∥x_{n+1}−x_{n}∥}}{2M}, \frac{1+\sqrt{1−4M∥x_{n+1}−x_{n}∥}}{2M}] .

Theorem. If we have 4M∥x_{n+1} − x_{n}∥ ≤ 1 and D contains the closed ball of radius

δ_{n} := [1−\sqrt{1−4M∥x_{n+1}−x_{n}∥}]/(2M) centered at x_{n}

then f admits a root x_{∗} ∈ D lying within the distance δ_{n} from x_{n}.

Corollary. Since Mδ_{n} ≤ 1/2 < 1 and ∥x_{n} − x_{∗}∥ ≤ δ_{n}, we have

∥x_{n+k} − x_{∗}∥ ≤ M^{2^{k}−1}δ_{n}^{2^{k}} ↘ 0 (k → ∞)

during the continuation of the Newton iteration.

CHAOS WITH NEWTON ITERATION

Recall. f : IR → IR C^{2}-smooth, f(x_{∗}) = 0, f′(x_{∗}) ≠ 0. Starting point x_{0}, x_{n+1}: location of the zero of the 1st-order Taylor polynomial of f around x_{n}

Heuristics: f ≈ [x ↦ f′(x_{n})(x − x_{n})] around x_{n}.

IN SEVERAL DIM as well: x_{n+1} := x_{n} − f′(x_{n})^{−1}f(x_{n}),

If f : IR^{N} → IR^{N} is C^{2}-smooth, the same notations can be used (FRÉCHET).

Discussion: Assume U is a convex neighborhood of x_{∗} and f′(x) is invertible for all x ∈ U.

Then

\begin{matrix}0 & = f(x_{∗}) \\ 0 & = f(x_{n}) + f′(x_{n})(xn + 1 − x_{n}) \end{matrix} \} \begin{matrix}f(x_{∗}) =^{h:=x_{∗}−x_{n}}= f(x_{n} + h) = \\ = f(x_{n})+f′(x_{n})h+∫_{t=0}^{1} ∫_{s=0}^{t}f′′(x_{n} + sh)h^{2} dsdt \end{matrix}

x_{n+1} − x_{∗} = f′(x_{n})^{−1} ∫_{t=0}^{1} ∫_{s=0}^{t}f′′(x_{n} + sh)(x_{n} − x_{∗})^{2} dsdt.

if x_{n} ∈ U and ,

∥x_{n+1} − x_{∗}∥ ≤ M_{U}∥x_{n} − x_{∗}∥^{2}, where M_{U} := \frac{1}{2} max_{x∈U} ∥f′(x)^{−1} max_{y∈U} ∥f′′(y).

Good starting guess: x_{0} ∈ M_{U} where M_{U}∥x_{0} − x_{∗}∥ < 1. Then induction yields that

∥x_{n} − x_{∗}∥ ≤ [M_{U}∥x_{0} − x_{∗}∥]^{2^{n}−1}∥x_{0} − x_{∗}∥ ↘ 0 (n = 1, 2, …).

Case of arbitrary starting point. We can encounter chaos even if looking for roots of polynomals (in 1 variable).

Observation. Regarding a holomorphic function f : C → C as a real mapping IR^{2} → IR^{2}, the steps of the Newton iteration can be written as

z_{n+1} := N(z_{n}), where N(z) := z − \frac{f(z)}{f′(z)} compl. diff.

Recall. (A celebrated theorem of Sharkowskii) The iterations of a non 3-periodic continuous self mappings of a connected compact topological space are chaotic provided it admits a fixed point of 3rd order. [Namely there is an invariant set being topologically equivalent to the Cantor set containing a dense family of fixed points of arbitrarily large orders.]

Let us apply the above mapping N(·) of the Newton iteration to a polynomial

f(z) := ∏_{k=1}^{N}(z − ω_{k}), ω_{k} ≠ ω_{ℓ} (k ≠ ℓ)

with 3 simple roots and extend it continuously to the Riemann sphere \overline{C} := C ∪ \{∞\} by setting N(∞), N(η_{1}), …, N(η_{N−1}) := ∞ at the roots η_{k} of the derivative N′(·). Given a fixed point z of 3rd order of N(·), we have

z = N^{∘3}(z) = N(N(N(z))), N(u) = u − \frac{1}{∑{N \atop k=1} \frac{1}{u−ω_{k}}} .

Since N(·) is a rational function, the function z ↦ z − N^{∘3}(z) is also rational. Thus we can write the equation for the fixed point of 3rd order in the polynomial form

0 = z − N^{∘3}(z) = \frac{P (z)}{Q(z)} ⟺ P (z) = 0.

This admits a complex solution by the Fundamental Theorem of the Algebra. Therefore N(·) admits a fixed point of 3rd order in any case. Thus we obtained:

Theorem. The Newton iteration algorithm applied to a complex polynomial of one variable with simple roots is either 3-periodic or chaotic with some starting point.

Remark. 3-periodicity can be checked with computer algebra if the coefficients are given numbers in the form algebraic formulas with integers. [There are not too many such cases].

The theoretical investigation of which cases are 3-periodic, is harder.

Example. Determine the unique real root of the polynomial f(x) := x^{3} − x − \frac{1}{\sqrt{3}} with Newton iteration. Thus

x_{n+1} = N(x_{n}), where N(x) := x − \frac{f(x)}{f′(x)} = x − \frac{x^{3} − x − \frac{1}{\sqrt{3}}}{3x^{2} − 1} (n = 0, 1, …).

Observation. Given any point x the half line (−∞, 0] × ×0 in the X-axis, we can draw exactly two tangent straight lines from x to the graph of f. Therefore

\begin{matrix}z < 0 & ⟹ N^{−1}\{z\} = \{x : N(x) = z\} = \{N_{−}^{−1}(z), N_{+}^{−1}(z)\}, \\ where N_{−}^{−1}(z) = [x < z : N(x) = z], N_{+}^{−1}(z) = [x > z : N(x) = z] . \end{matrix}

Since f is concave on the interval (−∞, 0] and its local maximum is attained at the point −1/\sqrt{3}, furthermore since the tangent straight line to the graph of f at the point of inflection (1/\sqrt{3}, f(1/\sqrt{3})) passes through the origin, we have

N_{−}^{−1} : IR ↔ IR, 0 ↦ −1/\sqrt{12} increasing, N_{+}^{−1} : (−∞, −1/\sqrt{3}] ↔ (−∞, 0] increasing and hence

I_{0} := [−1/\sqrt{3}, 0], I_{1} := N_{−}^{−1}1(I_{0}) ∩ (−∞, 0], I_{2} := N_{−}^{−1}(I_{1}) ∩ (−∞, −1/2] are well-defined intervals with the ordering I_{2} < I_{1} < I_{0}, furthermore N : I_{2} ↔ I_{1} ↔ I_{0} in increasing manner. The properties of the inverse N_{+}^{−1} on the right hand side give rise to the construction of the intervals

J_{0} := N_{+}^{−1}(I_{2}), J_{1} := N_{−}^{−1}(J_{0}), J_{2} := N_{−}^{−2}(J_{0}) = N_{−}^{−1}(J_{1})

where J_{k} ⊂ I_{k} (k = 1, 2, 3) and hence N : J_{0} ↔ J_{1} ↔ J_{2}. That is

J_{2} ⊂ I_{0}, N^{3} = N ∘ N ∘ N : J_{2} ↔ I_{2}.

Since the functions N, N^{2}, N^{3} are increasing on the interval I_{0} = [−1/\sqrt{3}, 0], the relation J_{2} ⊂ I_{2} implies that the function x − N^{3}(x) is ≤ 0 at the left end point of J_{2} and it is ≥ 0 at the right end point. Thus it admits a root in J_{2} by Bolzano’s classical theorem, that is exists a point z_{∗} ∈ J_{2}, which is a fixed point of 3rd order of the Newton iteration map N: (z_{∗} = N^{3}(z_{∗}) = N(N(N(z_{∗}))). This fact implies the (Sharkovskii type) chaotic behaviour of the Newton iteration.

num0 eng fig1.pdf
3 pages

Figure by Prof. R. Vajda Figure by Prof. R. Vajda constructed with Wolfram Mathemaica 12 Figure by Prof. R. Vajda constructed with Wolfram Mathemaica 12 Figure by Prof. R. Vajda constructed with Wolfram Mathemaica 12 Figure by Prof. R. Vajda constructed with Wolfram Mathemaica 12 Figure by Prof. R. Vajda constructed with Wolfram Mathemaica 12

CONDITION NUMBERS

Recall. The number x̃ ∈ IR approaches the number x ∈ IR with (absolute) error |x̃ − x|, with relative error |x̃ − x|/|x|. In general: In a normed space (X, ∥.∥), the point x̃ ∈ X approaches the point x(∈ X) with absolute resp. relative error]

err_{x}(x̃) := ∥x̃ − x∥, rel_{x}(x̃) := \frac{∥x̃ − x∥}{∥x∥}.

Definition. Let U open ⊂ X, and f : U → X. The condition number of the function (mapping) f at the point x ∈ U resp. on the region U ⊂ X is

cond_{x}(f) := lim sup_{_{x}→x} \frac{rel_{f(x)}(f(x̃))}{rel_{x}(x̃)}, cond_{U}(f) := sup_{x∈U} cond_{x}(f).

Lemma. If f is C^{1}-smooth then cond_{x} = ∥x∥ · ∥f′(x)∥/∥f(x)∥.

Proof. cond_{x}(f) = lim sup_{v→0} \frac{∥f(x + v) − f(x)∥/∥f(x)∥}{∥(x + v) − x∥/∥x∥} =

= lim sup_{v→0} = \frac{∥f′(x)v + o(∥v∥)∥/∥f(x)∥}{∥v∥/∥x∥} =

= lim sup_{v→0} \frac{∥f′(x)v∥}{∥v∥} · \frac{∥x∥}{∥f(x)∥} = \frac{∥f′(x)∥ · ∥x∥}{∥f(x)∥}.

Condition number of a matrix. As usually, we identify the real (N × N)-es matrix with the linear mapping x ↦ (IR^{N}→IR^{N}). In such manner, if x ≠ 0 and 0 ≠ U then

cond_{x}(A) = lim sup_{v→0} \frac{∥Av∥}{∥v∥} · \frac{∥x∥}{∥Ax∥} = ∥A∥\frac{∥x∥}{∥Ax∥}.

Thus if A is invertible and U is an arbitrary neighborhood of 0 then

cond_{U\\{0\}}(A) = sup_{x≠0} cond_{x}(A) = ∥A∥ sup_{x≠0} \frac{∥x∥}{∥Ax∥} =^{y=Ax}

= ∥A∥ sup_{y≠0} \frac{∥A^{−1}y∥}{∥y∥} = ∥A∥ · ∥A^{−1}∥.

Definition. The condition number of the matrix A is

cond(A) := ∥A∥ · ∥A^{−1}∥ .

We obtain with direct calculation:

Proposition. (1) cond_{x}(f + g) ≤ cond_{x}(f) + cond_{x}(g),

(2) cond_{x}(f ∘ g) ≤ cond_{g(x)}(f) · cond_{x}(g),

(3) 1 ≥ cond(A) = cond(A^{−1}) ≤ cond(B) · cond(C) ha A = BC,

(4) With respect to L^{2}-norm, if Q is an orthogonal matrix then cond(Q) = 1 and cond(QX) = cond(XQ) = cond(X).

LINEAR ERROR ANALYSIS

Basic setting: A, δA ∈ Mat(N, N, IR), b, δb, x, δx ∈ Mat(N, 1, IR) are given matrices resp. column vectors (here δ is no operator, but only a notation for the ”error term”) where

Ax = b (is the ideal equation),

(A + δA)(x + δx) = b + δb (is the computed solution).

Assumption (technical): det(A) ≠ 0, ∥A^{−1}δA∥ < 1.

Remark. In the basic setting, A + δA is invertible. Namely

(A + δA)^{−1} = [A(1 + A^{−1}δA)]^{−1} = (1 + A^{−1}δA)^{−1}A^{−1} =

= ∑_{n=0}^{∞}(−1)^{n}[A^{−1}δA]^{n}A^{−1},

∥(A + δA)^{−1}∥ ≤ ∑_{n=0}^{∞} ∥A^{−1}δA∥^{n}∥A^{−1}∥ = \frac{∥A^{−1}∥}{1 − ∥A^{−1}δA∥}.

Recall: cond(A) = ∥A∥ ∥A^{−1}∥ is the condition number of the matrix A.

Fundamental Theorem (on data sensibility). Under the above hypothesis, we can estimate the relative errors as

\frac{∥δx∥}{∥x∥} ≤ \frac{cond(A)}{1 − ∥A^{−1}δA∥} (\frac{∥δA∥}{∥A∥} + \frac{∥δb∥}{∥b∥}) .

Proof. (A + δA)(x + δx) = b + δb, Ax = b, ⟹ \underline{(δA)x + (A + δA)δx = δb},

δx = (A + δA)^{−1}[δb − (δA)x],

∥δx∥ ≤ ∥(A + δA)^{−1}∥[∥δb∥ + ∥δA∥ ∥x∥],

\frac{∥δx∥}{∥x∥} ≤ ∥(A + δA)^{−1}∥ [\frac{∥δb∥}{∥x∥} + ∥δA∥] .

Observation: \frac{1}{∥x∥} = \frac{1}{∥A^{−1}b∥} ≤ \frac{1}{∥b∥ inf\{∥A^{−1}e∥ : ∥e∥ = 1\}} =

= \frac{1}{∥b∥} \frac{1}{inf\{∥A^{−1}y∥/∥y∥ : y ≠ 0\}} =

= \frac{1}{∥b∥} sup\{∥y∥/∥A^{−1}y∥/ : y ≠ 0\} =

= \frac{1}{∥b∥} sup\{∥Az∥/∥z∥/ : z ≠ 0\} = \frac{1}{∥b∥}∥A∥.

Thus

\begin{aligned}\frac{∥δx∥}{∥x∥} & ≤ ∥(A + δA)^{−1}∥ (\frac{∥δb∥}{∥b∥}∥A∥ + ∥δA∥) ≤ \\ & ≤ \frac{∥A^{−1}∥}{1 − ∥A^{−1}δA∥}∥A∥ (\frac{∥δA∥}{∥A∥} + \frac{∥δb∥}{∥b∥}) . Qu.e.d. \end{aligned}

Corollary. ∥A^{−1}∥ ∥δA∥ < 1 ⟹

\frac{∥δx∥}{∥x∥} ≤ \frac{cond(A)}{1 − cond(A)∥δA∥/∥A∥} (\frac{∥δA∥}{∥A∥} + \frac{∥δb∥}{∥b∥}) .

Remark. In our above considerations, we can use any matrix norm defined by any (vector) norm ∥ · ∥ on IR^{N} ≡ Mat(N, 1, IR).

We apply mostly the Euclidean norm

∥z∥_{2} := [ ∑_{k=1}^{N} z_{k}^{z}]^{1/2} (for vectors),

∥B∥_{2} := sup {∥Bz∥_{2} : ∥z∥_{2} = 1} (for matrices).

It is an interesting geometrical fact that

\frac{1}{cond^{∥·∥_{2}}(A)} = min \{\frac{∥δA∥_{2}}{∥A∥_{2}} : A + δA is not invertible\} .

Lower estimate for \frac{∥δx∥}{∥x∥}: δA = 0 ⟹ \frac{1}{cond(A)} \frac{∥δb∥}{∥b∥} ≤ \frac{∥δx∥}{∥x∥}.

Proof. δA = 0 ⇒ Aδx = δb ⇒ ∥δx∥ ≥ ∥δb∥/∥A∥. Since ∥x∥ = ∥A^{−1}b∥ ≤ ∥A^{−1}∥ ∥b∥, we obtain

∥δx∥/∥x∥ ≥ [∥δb∥/∥A∥]/[∥A^{−1}∥ ∥b∥] = ∥δb∥/[cond(A)∥b∥].

Estimation componentwise in case of (δb = 0 resp. δA = 0). Notation. |Z| := [|z_{ij}|]_{i=1}^{K}_{j=1}^{M} for any K × M matrix; e^{[i]} := [i. unit vector].

Theorem. Assume we have

|δA| ≤ γ|A|, |δb| ≤ γ|b|

in the equation (A + δA)(x + δx) = b + δb on data sensibility. Then, in terms of the column vectors

s^{[i]} := [e^{[i]}]^{T}A^{−1}

of the matrix A^{−1} we can conclude the following estimates.

(1) In case of δb = 0 and given x̂ := x + δx we have |δx_{i}| ≤ γ|s^{[i]}| |A| |x̂|;

(2) In case of δA = 0 we have |δx_{i}| ≤ γ|s^{[i]}||b|, \frac{|δx_{i}|}{|x_{i}|} ≤ γ \frac{|s^{[i]}||b|}{|s^{[i]}b|}.

Proof. (A + δA)(x + δx) = b + δb ⇒

Aδx = −(δA)(x + δx) + δb,
δx = A^{−1}[−(δA)(x + δx) + δb],
δx_{i} = [e^{[i]}]^{T}δx = [e^{[i]}]^{T}A^{−1}[−(δA)(x + δx) + δb] =
= s^{[i]}[−(δA)(x + δx) + δb].

(1) IF δb = 0, |δx_{i}| ≤ |s^{[i]}| |δA| | \underbrace{x + δx}_{_{x}} | ≤ γ|s^{[i]}| |A| |x̂|;

(2) IF δA = 0, |δx_{i}| ≤ |s^{[i]}| |δb| ≤ γ|s^{[i]}| |b|,

\frac{|δx_{i}|}{|x_{i}|} ≤ γ \frac{|s^{[i]}||b|}{|e^{[i]}A^{−1}b|} ≤ γ \frac{|s^{[i]}||b|}{|s^{[i]}b|}. Q.e.d.

A priori analysis backwords Basic setting: x̂=x+δx=Cb. E.g. C =[A^{−1} with rounding errors] not computed explicitly; A priori assumption: C = (A + δA)^{−1}, where δA not know for us. Lemma. If the norm of R := I − AC is < 1 then C is invertible and

∥A^{−1}∥ ≤ \frac{∥C∥}{1 − ∥R∥}, \frac{∥R∥}{∥C∥} ≤ ∥C − A^{−1}∥ ≤ \frac{∥C∥ ∥R∥}{1 − ∥R∥}.

Proof. R = I − AC ⇒ C = A^{−1}(I − R) is a product of invertible matrice if, ∥R∥ < 1.

R=I −AC ⇒ A=(I −R)C^{−1} ⇒ A^{−1} =C(I −R)^{−1} =C ∑_{n=0}^{∞}R^{n} ⇒ ∥A^{−1}∥≤∥C∥ ∑_{n=0}^{∞} ∥R∥^{n}.

R=I −AC ⇒ R=A(A^{−1}−C) ⇒ ∥R∥≤∥A∥ ∥A^{−1}−C∥.

Here we have A^{−1}−C =C(I −R)^{−1}−C = C ∑_{n=0}^{∞}R^{n} ⇒ ∥A^{−1}−C∥≤∥C∥/(1 − ∥R∥).

Lemma. In case of C = Â^{−1} = (A + δA)^{−1} and ∥R∥ < 1 we have

∥δA∥ ≤ ∥R∥ ∥C^{−1}∥ ≤ \frac{∥R∥ ∥A∥}{1 − ∥R∥}.

Proof. Observation: δA = C^{−1} − A = (I − AC)C^{−1} = RC^{−1}. Here we have C^{−1} = [A^{−1}(I − R)]^{−1} = (1 − R)^{−1}A whence the norm estimates are immediate. A posteriori analysis Basic setting: Given an approximate solution

y(= x + δx) of the equation Ax = b (with det(A) ≠ 0).

Let us estimate the residual- resp error vectors

r := Ay − b( = (δA)(x + δx) − δb), e := A^{−1}r.

On the basis of our earlier results we can conclude the following.

Lemma. ∥e∥ ≤ ∥r∥ ∥A^{−1}∥ ≤ \frac{∥r∥ ∥C∥}{1 − ∥R∥}, \frac{∥e∥}{∥x∥} ≤ cond(A)\frac{∥r∥}{∥b∥}.

On the effect of the inaccuracy in machine representations Let r̂ = r + δr := [machine output for r(= b − Ay)]. Assumption. A ∈ IR^{N×N}, and by writing u for the machine-0,

|δr| ≤ γ_{N+1}(|A| |y| + |b|), ahol γ_{N+1} := \frac{(N + 1)u}{1 − (N + 1)u}.

Theorem. Then with the norm ∥.∥ := ∥.∥_{∞} we have

\frac{∥e∥_{∞}}{∥y∥_{∞}} ≤ \frac{∥|A^{−1}|(|r̂| + γ_{n+1}(|A| |y| + |b|))∥_{∞}}{∥y∥_{∞}} .

The effect of rounding in LU-factorization Machine rounding: The machine output of the operations ♢(= +, −, ·, /) is

[x♢y] = (1 + δ_{♢,x,y})(x♢y), |δ_{♢,x,y}| ≤ u := [machine 0].

Assumption. (A + δA) = (L + δL)(U + δU) where A = LU, L is a lower-triangular matrix with main diagonal 1, U is an upper-triangular matrix, both of N × N type, the machine outputs are Â:= A + δA, L̂ := L + δL, Û := U + δU where u := [machine-0].

One can establish the following fact. Theorem. Under the above assumption,

A = LU, Â = L̂Û, [main diagonal of L̂ ] = 1 ⟹ |δA| ≤ \frac{Nu}{1 − Nu}|L̂| |Û|.

Excercise: verify the numerical solution for a system of linear equations Theoretical problem: Ax = b. We are given a numerical solution x̃. That is

Ãx̃ = b̃, x = x̃ + e, A = Ã + E, b = b̃ + d,

where we have estimates for the errors d resp. E. In the most general setting

d ∈ D, E ∈ E,

where the family D of vectors and the family E of matrices is given. Question: How much is x̃ suitable instead of x? Estimate in terms of D, E.

(Ã + E)(x̃ + e) = b̃ + d,
Ex̃ + Ee + Ãe = d,
e = (Ã + E)^{−1}(d − Ex̃).

Starting argument: e ∈ {(Ã + E)^{−1}(d − Ex̃) : E ∈ E, d ∈ D}. Example. D := {d′ : ∥b′∥ < δ}, E := {E′ : ∥E′∥ < ε}. Ekkor

∥e∥ ≤ sup_{∥E∥<ε} ∥(Ã + E)^{−1}∥(δ + ε∥x̃∥).

In general ∥B^{−1}∥ = sup_{∥z∥=1} ∥B^{−1}z∥ = \frac{1}{inf_{∥y∥=1} ∥By∥}. Hence

∥(Ã + E)^{−1}∥ ≤ \frac{1}{[inf_{∥y∥=1} ∥Ãy∥ − ε]_{+}} = \frac{1}{\frac{1}{∥_{A^{−1}}∥} − ε},

provided ∥Ãy∥ > ε (∥y∥ = 1). This latter can be rewritten in terms of te conditional number

K(B) := ∥B^{−1}∥ ∥B∥

whence

\begin{matrix}(∥(Ã + E)^{−1}∥ ≤ \frac{Ã^{−1}}{1 − ε∥(Ã + E)^{−1}∥} = \frac{K(Ã)}{∥Ã∥ − εK(Ã)}, \\ ∥e∥ ≤ \frac{K(Ã)}{∥Ã∥ − εK(Ã)} (δ + ε∥x̃∥). \end{matrix}

f ≺ F: ∀ x ∈ dom(f) ∃ I ∈ dom(F) and f(x) ∈ F (I) whenever x ∈ dom(f) with x ∈ I ∈ dom(F). F_{1} ≺ F_{2}: dom(F_{1}) ⊂ dom(F_{2}) and F_{1}(I) ⊂ F_{2}(I) whenever I ∈ dom(F_{1}) dist(x, Y) = inf{d(x, y) : y ∈ Y } dist(X, Y) = [Hausdorff distance] = sup {d(x, Y), d(y, X) : x ∈ X, y ∈ Y } = = max { sup{d(x, Y) : y ∈ Y \ X}, sup{d(y, X) : x ∈ X \ Y } If I = [a, b], J7[c, d] finite intervals, then

(1) a ≤ c ≤ b ⇒ dist(I, J) = max {|c − a|, |b − d|},

(2) b ≤ c ⇒ dist(I, J) = max {|a − d|, |b − c|}. diam(X) = sup{d(x, y) : x, y ∈ X} Lemma. f ≺ F, I ⊂ I_{1} ∪ ⋯ ∪ I_{M}, ⇒ f(I) ⊂ ⋃_{k=1}^{M} F (I_{k}), dist(range(f|I), ⋃_{k=1}^{M} F (I_{k})) ≤ max_{k=1}^{M} diam(F (I_{k})) Lip(f) = sup_{x,y∈dom(f)} \frac{|f(x)−f(y)|}{|x−y|} Lip(F) = sup_{I∈dom(F)} \frac{diam F (I)]}{diam I} X ⊂ dom(f) ⇒ diam f(X) ≤ Lip(f) · diam X I ∈ dom(F) ⇒ diam F (I) ≤ Lip(F) · diam I f ∈ C^{1}(I) ⇒ Lip(f) = sup |f′(I)| Definition. An interval function F : X → Intv(IR) is isotonic (F ↗ in notation) if I ∈ X ⇒ ∀ J ∈ Intv(I) J ∈ X and F (J) ⊂ F (I) Lemma. F_{1}, F_{2} isotonic ⇒ F_{1} ∘ F_{2} isotonic. Evaluation. Φ ⟶ x = Φ_{0}, Φ_{1}, …, Φ_{n} = Φ P hi_{i} = [x OR const. OR s(Φ_{j}) for some j < i, s ∈ 𝔖 OR Φ_{j}♢P hi_{k} for some j, k < i and operation ♢] fΦ_{0} = x, ^{f}P hi_{i} = [x RESP. const. RESP. s(^{f}Φ_{j}) RESP. Φ_{j}♢P hi_{k}] 𝔖Φ_{0} = X = [ξ ↦ ξ], 𝔖P hi_{i} = [X RESP. Const. RESP. S(𝔖Φ_{j}) RESP. 𝔖Φ_{j}♢𝔖P hi_{k}]

INTERVAL ARITHMETICS

Basic concepts. As usually, by real intervals we mean the subsets (a, b), [a, b), (a, b], [a, b] with a ≤ b ∈ IR ∪ {±∞}. The family of all subintervals of the interval I will be denoted by Intv(I) := {J interval: J ⊂ I}. In order to control the accuracy of calculating the value of a formula Φ at a given argument x, like Φ(x) = \frac{e^{x}}{sin^{2}(x) + 2 cos^{2} x}, we consider x as belonging to some (open) interval I containing the set of all possible values of the representations of x with error (like 3.14 fo π) and, for each elementary step of the calculation, try to find an interval which inclues the partial result of the step. In general, if X_{1}, …, X_{n}, Z are arbitrary (non-empty) sets and f : X_{1} × ⋯ × X_{n} → Z is a mapping of n variables, we extend the effect of f to sets X_{k} ⊂ X_{k} in the natural manner as f(X_{1}, …, X_{n}) := {f(x_{1}, …, x_{n}) : x_{k} ∈ X_{k} (k = 1, …, n)}. The classical algebraic operations ♢ = +, −, ·, / are regarded as functions of two variables: X♢Y := {x♢y ∈ IR ∪ {±∞} : x ∈ X, y ∈ Y } (X, Y ⊂ IR ∪ {±∞}). Inclusion Principle. X_{1} ⊂ Y_{1}, …, X_{n} ⊂ Y_{n} ⟹ f(X_{1}, …, X_{n}) ⊂ f(Y_{1}, …, Y_{n}). Theorem. Intv♢Intv = Intv (♢ = +, −, ·, / (0 ∉ denominator)). In any case we can write [a, b]♢[c, d] in the form [a, b]♢[c, d] = [φ(a, b, c, d), ψ(a, b, c, d)]. Proof. Recall: continuous functions map compact connected sets into compact connected sets and the compact connected subsets of IR are exactly the closed bounded intervals. Also the operations ♢ = +, −, ·, / (0 ∉ denominator) are continuous and monotonic in their variables over the plane quartals IR_{ε} × IR_{σ} (ε, σ = ±). Hence we can write the end points as formulas of the end points a, b, c, d of the two intervals.

Example. [a, b]♢[c, d] = [min{a♢c, a♢d}, max{b♢c, b♢d}] (♢ = ±). Exercise. 1) End points of [a, b] · [c, d]. 2) End points of [a, b]/[c, d].

Remark. We can define a (bit sophisticated) arc arithmetics if we identify the unit circle with the one point compactification IR∪{∞} of IR, the compact connected subsets are identitified with the closed arcs. Advantage: continuous maps send intervals to intervals without any restriction in this setting. Nevertheless we shall not discuss this model in our course.

In the sequel we need extensions for the inclusion property I ⊂ J ⇒ f(I) ⊂ f(J). In general, let X, Z ≠ ∅ arbitrary sets and F : X → Z be a set mapping where X ⊂ 2^{X}(:= {subsets ofX}) and Z ⊂ 2^{Z}. We say that F is isotonic if

F ↗: X, Y ∈ X and Y ⊂ X ⟹ F (Y) ⊂ F (X). Lemma. F_{1}, F_{2} isotonic and rangeF_{1} ⊂ dom(F_{2}), ⇒ F_{1} ∘ F_{2} isotonic.

By an interval function we mean a map F : X → Intv[−∞, ∞] where X ⊂ 2^{[−∞,∞]} with hereditary domain i.e. such that Intv(J) ⊂ X whenever J ∈ X.

Lemma. The family of isotonic interval functions is closed for composition. (Notice: dom(F_{1} ∘ F_{2}) = {I : F_{2}(I) ∈ dom(F_{1})}.)

In general, an ordinary function f : X ⊃ dom(f) → Z is included in the set function F : 2^{X} ⊃ dom(F) → 2^{Z} if

f ≺ F: ∀ x ∈ dom(f) ∃ X ∈ dom(F) f(x) ∈ F (X). Given two set functions F_{k} : 2^{X} ⊃ dom(F_{k}) → 2^{Z} k = 1, 2), we also introduce the concept of inclusion as

F_{1} ≺ F_{2}: dom(F_{1}) ⊂ dom(F_{2}) and F_{1}(I) ⊂ F_{2}(I) whenever I ∈ dom(F_{1}) Remark. We intend to use inclusions of a function by a complicated formula with interval functions of simpler formulas on small intervals such that the graphs do not differ too much. This is the trick of classical mathematics in several cases when establishing estimates by constants or simple formulas. With computer, we can enhance this old techniques enormously.

Recall. In a metric space (X, d), the diameter of a set A ⊂ X resp. the distance of a point x ∈ X from a set A ⊂ X are

dist(x, A) := inf{d(x, a) : a ∈ A.}. We can use the Hausdorff distance dist(A, B) for measuring the difference between the sets A, B ⊂ X wrt. the underlying distance:

dist(X, Y) = [Hausdorff distance] = sup {d(x, Y), d(y, X) : x ∈ X, y ∈ Y } =

= max { sup{d(x, Y) : y ∈ Y \ X}, sup{d(y, X) : x ∈ X \ Y }.

Excercise. If I = [a, b], J = [c, d] finite intervals, then

(i) a ≤ c ≤ b ⇒ dist(I, J) = max {|c − a|, |b − d|}, (ii) b ≤ c ⇒ dist(I, J) = max {|a − d|, |b − c|}.

Lemma. Covering Lemma). Assume f ≺ F and X ⊂ X_{1} ∪ ⋯ ∪ X_{M} (where X ∈ dom(f) and X_{k} ∈ dom(F)). Then (1) f(X) ⊂ F (X_{1}) ∪ ⋯ ∪ F (X_{n}),

(2) dist(range(f|X), ⋃_{k=1}^{M} F (X_{k})) ≤ max_{k=1}^{M} diam(F (X_{k})) .

Proof. (1) Given x∈X, x∈X_{k} for some k =k(x). Hence f(x)∈F (X_{k(x)})⊂⋃_{m=1}^{M} F (X_{m}).

(2) Since, by (1), f(X) ⊂ ⋃_{m=1}^{M} Y_{m}, we simply have dist(f(X), ⋃_{m=1}^{M} Y_{m}) = sup {dist(f(x), ⋃_{m=1}^{M} Y_{m}) : x ∈ X}.

Let us write Y_{m} := F (X_{m}) with f(x) ∈ Y_{k(x)} (x ∈ X) and for short. Observe that

dist(f(x), ⋃_{m=1}^{M} Y_{m}) ≤ dist(f(x), Y_{k(x)}) (x ∈ X) ≤ diam(Y_{k(x)}. Hence the statement is immediate.

Algebraic formulas with interval arithmetics

Henceforth we fix a family

S := {x, sin, cos, log, … }

of classicaal real (elementary) functions. We shall call the members of S special functions. We also fix an isotonic interval extension to each special function for the family

𝔖 := {X, Sin, Cos, Log, … }

Remark. The family 𝔖 will play the role of the machine representation of the interval functions I ↦ I, I ↦ sin(I) = {sin x : x ∈ I} etc. Classical computer arithmetics work only with finitely many number symbols (e.g. with flooting point numbers of 48 digits in many practical cases). In such machines the intervals X(I), Sin(I) etc. are represented only with endpoints from the given number symbols. Therefore we do not assume automatically that X(I) = I, Sin(I) = {sin x : x ∈ I} etc.. Of course X(I) ⊂ I, Sin(I) ⊂ {sin x : x ∈ I} etc. in any case.

Definition. Let Φ = Φ(x) denote an elementary function expression of the variable symbol x. That is Φ consists of finitely many function symbols from S combined with the operations ♢ ∈ {±, ·, /, max, … } in a syntactically correct manner. (Example: Φ = Φ(x) = log(log x)/\sqrt{2 + sin^{3} x}).

Evaluation. The usual mathematical form is a shorthand way of describing the process how Φ is evaluated. We identify Φ with a sequence

x = Φ_{0}, Φ_{1}, Φ_{2}, …, Φ_{n} = Φ such that

Φ_{i} =[x OR const. OR s(Φ_{j}) for some j <i, s∈𝔖 OR Φ_{j}♢P hi_{k} for some j, k < i]. We introduce the natural numeric- resp. interval realizations^{f}Φ resp. ^{I}Φ of Φ with respect to the representation 𝔖 of the special functions as the numeric- resp. interval functions

fΦ(x) := [Formal step by step substitution of x into Φ],

𝔖Φ(J) := [Formal step by step substitution wrt. 𝔖 of J into Φ].

That is, in terms of the evaluation sequence, fΦ_{0} = x, ^{f}Φ_{i} = [x RESP. const. RESP. s(^{f}Φ_{j}) RESP. Φ_{j}♢Φ_{k}]

𝔖Φ_{0} =[ξ ↦ ξ], 𝔖Φ_{i} =[X RESP. Const. RESP. S(𝔖Φ_{j}) RESP. 𝔖Φ_{j}♢𝔖Φ_{k}].

Example. 1) For Φ = log(log x)/\sqrt{2 + sin^{3} x}, Φ_{0} = x, Φ_{1} = log(x), Φ_{2} = log(log(x)), Φ_{3} = 2, Φ_{4} = x, Φ_{5} = sin(x), Φ_{6} = (sin(x))^{3}, Φ_{7} = Φ_{3} + Φ_{6}, Φ_{8} = \sqrt{Φ_{7}}, Φ_{9} = Φ_{2}/Φ_{8} = Φ.

2) For Ψ_{1}(x) := 1 − x · x, Ψ_{2}(x) := (1 − x) · (1 + x) we have ^{f}Ψ_{1} = ^{f}Ψ_{2} : x ↦ 1 − x^{2}

\begin{matrix}𝔖Ψ_{1}[−1, 1] = 1 − [−1, 1] · [−1, 1] = [1, 1] − [−1,] = [0, 2], \\ 𝔖Ψ_{2}[−1, 1] = (([1, 1] − [−1, 1]) · (([1, 1] + [−1, 1]) = [0, 2] · [0, 2] = [0, 4]. \end{matrix}