Numerical analysis - Appunti + esami svolti

Numerical analysis notes

Lorenzo Sostegni

Contents

• Foundations of numerical mathematics
• Polynomial approximation of functions and data
• Numerical linear algebra
• Ordinary differential equations
• Nonlinear equations

Foundations of numerical mathematics

1. Definition of mathematical problem

Let D be the set of input data for the model and let V be the set of solutions. Let D and V be vector spaces. The mathematical problem reads: given d ∈ D, find x ∈ V such that F(x, d) = 0. The mapping F defined as F: V × D → Y, (x, d) ↦ F(x, d), is referred to as the mathematical problem.

2. Well-posedness of the mathematical problem

The mathematical problem is said to be well-posed if:

• Existence: given any d ∈ D there exists x ∈ V such that F(x, d) = 0
• Uniqueness: given any d ∈ D, if there exist x, y ∈ V such that F(x, d) = 0 and F(y, d) = 0, then x = y
• Continuous dependence on data: given any d, d + δd and their corresponding unique solutions x, x + δx ∈ V such that F(x, d) = 0 and F(x + δx, d + δd) = 0, there exist two positive constants η and K, which may depend on d such that if ||δd|| ≤ K(d) ||δx|| ≤ η(d)||δd||, then d ∈ D.

3. Stability of the mathematical problem

The property of continuous dependence on data is also referred to as the property of stability. Indeed, the fact that the solution depends continuously on the data means that a small perturbation δd in the input data d for the model leads to a small perturbation δx in the solution x. In other words, if we change the input data only slightly, the solution will not change that much, meaning that the problem is stable.

4. Definition of numerical problem

Let D_h and V_h be finite-dimensional subspaces of D and V, respectively. The numerical problem reads: given d_h ∈ D_h, find x_h ∈ V_h such that F_h(x_h, d_h) = 0. The mapping F_h defined as F_h: V_h × D_h → Y, (x_h, d_h) ↦ F_h(x_h, d_h), is referred to as the numerical problem.

5. Well-posedness of the numerical problem

The numerical problem is said to be well-posed if:

• Existence: given any d_h ∈ D_h, there exists x_h ∈ V_h such that F_h(x_h, d_h) = 0
• Uniqueness: given any d_h ∈ D_h, if there exist x_h, y_h ∈ V_h such that F_h(x_h, d_h) = 0 and F_h(y_h, d_h) = 0, then x_h = y_h
• Continuous dependence on data: given any d_h, d_h + δd_h and their corresponding unique solutions x_h, x_h + δx_h ∈ V_h such that F_h(x_h, d_h) = 0 and F_h(x_h + δx_h, d_h + δd_h) = 0, there exist two positive constants η_h and K_h, which may depend on d_h such that if ||δd_h|| ≤ K_h(d_h) ||δx_h|| ≤ η_h(d_h)||δd_h||, then d_h ∈ D_h.

6. Residual

The residual associated with the numerical model is defined as: R(x, d) := F(x, d) - F_h(x, d).

7. Consistency of the numerical model

The numerical model is said to be consistent if lim_h→0 d_h = d or equivalently lim_h→0 R_h(x, d) = 0. Moreover, the numerical model is said to be consistent of order p with respect to the discretization parameter h if there exists C > 0 independent of h and p > 0 such that ||R_h(x, d)|| ≤ C h^p.

8. Definition of the discretization error

The error e_h := x - x_h is called the discretization error. Assume that there exists a positive constant C independent of h and a positive real number p such that ||e_h|| ≤ C h^p. In such a case, we say that x_h converges to x with order p with respect to h. Using a bilogarithmic scale, we can see that log(||e_h||)_V = log(C) + p log(h): the slope of the line gives us the order of convergence.

9. Absolute and relative errors

Absolute error: E_abs(x_h) := ||x - x_h||_V Relative error: E_rel(x_h) := ||x - x_h||_V / ||x||_V

10. Equivalence theorem (Lax-Richtmyer)

The numerical problem is convergent with respect to the discretization parameter h if and only if it is stable and consistent. Stability + Consistency ⇔ Convergence.

11. Machine epsilon

The machine epsilon is the smallest floating point number such that 1 + ε_M > 1. In our computers β = 2 and t = 53, so we have: ε_M = β^(1-t) = 2^-53 ≈ 2.2204 × 10^-16.

12. Forward, backward and centered finite differences

Forward difference: y'(t̄) ≈ (y(t̄ + h) - y(t̄)) / h Backward difference: y'(t̄) ≈ (y(t̄) - y(t̄ - h)) / h Centered difference: y'(t̄) ≈ (y(t̄ + h) - y(t̄ - h)) / (2h)

Polynomial approximation of functions and data

1. Global polynomial interpolation

In order to interpolate our function f(x): [a, b] → R, we need to select a polynomial y(x) ∈ P_n. The domain is divided in M subintervals, individuated by n + 1 distinct nodes. We need to satisfy:

    a₀ + a₁x₀ + a₂x₀² + ... + anx₀n = y₀
    a₀ + a₁x₁ + a₂x₁² + ... + anx₁n = y₁
    ...
    a₀ + a₁xₙ + a₂xₙ² + ... + anxₙn = yₙ

This system can be also seen as:

    X · a = y

Where the matrix X is called Vandermonde matrix, the vector a is the unknown datum and the vector y is the given datum. In order to solve this linear system, we need the Vandermonde matrix to be non-singular, and this always holds since we have n + 1 distinct nodes. We can collect the interpolating polynomial coefficients by simply computing a = X^-1 y.

2. A priori interpolation error estimate

For any x ∈ [a, b] the interpolation error is defined as E_nf(x) := f(x) - Π_nf(x).

||E_nf||_C([a,b]) ≤ (||f⁽ⁿ⁺¹⁾||_C([a,b]) / (n + 1)!) * ||ω_n+1||_C([a,b])

Where ω_n+1(x) = (x-x₀)(x-x₁)...(x-x_n) is the nodal polynomial.

3. Convergence of the interpolation

We say that Π_nf converges uniformly if lim_n→∞ Π_nf(x) = f(x) ∀x ∈ [a, b], in other words, lim_n→∞ E_nf(x) = 0.

lim_n→∞ ||E_nf||_C([a,b]) = 0 if ω_n+1(x) = 0 for all x ∈ [a, b].

4. Chebyshev's nodes

In this particular case, the domain is divided into subdomains of different lengths and in particular every node abscissa is expressed as:

x(x̄) = A x̄ + B cos(πi/n)

In this way, we have very tiny intervals at the edges of the domain and we reduce the interpolation error.

5. Piecewise polynomial interpolation

To construct the method, we introduce a partition τ_h of [a, b] into M intervals. We denote by h_i = x_i+1 - x_i the length of each interval and set h = max(h_i). We associate with τ_h the following space of piecewise polynomials:

X_h^r(τ_h) := {v ∈ C([a,b]) | v|_Ki ∈ P_r(K_i) ∀K_i ∈ τ_h}

dim(X_h^r) = rM + 1 is usually referred to as the finite element space of degree r associated with τ_h.

6. A priori piecewise interpolation error estimate

Let us assume that f ∈ C^r+1([a,b]), r ≥ 1. Then there exists a positive constant C independent of h such that:

||E_hf||_C([a,b]) ≤ C ||f^(r+1)||_C([a,b]) h^r+1

7. Lagrange polynomials

Lagrange polynomials are in the form:

L_i(x) = Π_{j=0, j≠i}^n (x - x_j) / (x_i - x_j)

Using these particular polynomials, we can interpolate our function with the following formula:

Π_n(x) = Σ_i=0^n y_i L_i(x)

8. Least-squares approximation

We want to find a polynomial p(x) such that:

Σ_i=0^m (y_i - p(x_i))² = min

The polynomial is found in the following way for p(x) = c₀ + c₁x:

    c0 + c1x0 = y0
    c0 + c1x1 = y1
    ...
    c0 + c1xm = ym

In algebraic terms:

    XTX · c = XTy

To solve this system, we multiply both terms by X^T and we obtain c = (X^TX)^-1X^Ty.

9. Quadrature rules

The error we make with quadrature rule is equal to:

    |E(f)| ≤ Ck hr+1 ||f(r+1)||C([a,b])

Midpoint Quadrature Rule: Q_i0,k(f) = f(x̄_i)h E(f) = (h³/24)f''(ξ_i)

Trapezoidal Quadrature Rule: Q_i1,k(f) = [f(x_i) + f(x_i+1)]h/2 E(f) = -(h³/12)f''(η_i)

Gauss-Legendre Quadrature Rule: x_GL,1 = x̄_i - √(3/2)h, x_GL,2 = x̄_i + √(3/2)h Q_GL,1,k(f) = (f(x_GL,1) + f(x_GL,2))h/2 |E(f)| = C₍₄₎(x) h⁵ ||f||_C(k)

Cavalieri-Simpson Quadrature Rule: Q_i2,k(f) = (f(x_i) + 4f(x̄_i) + f(x_i+1))h/6