1 About Error link

2 Root Finding link

2.1 Bisection Method link

Necessary Condition: $f$ is a continuous function $f \in C[a,b]$ and $f(a)f(b) < 0$.

2.1.1 Algorithms link

Given $f(x)$, $f(a)$, $f(b)$ with tolerance $\epsilon_1$, $\epsilon_2$:

def bisection(f, a, b, tolerance):
    if f(a) * f(b) > 0:
        return None, None

    if a > b:
        a, b = b, a

    mid = (a+b)/2
    iteration = 0
    error = [(b-a)/2]

    while (abs(f(mid)) > tolerance and b - a > tolerance):
        mid = (a + b)/2

        if f(mid) == 0:
            return mid, None
        elif f(mid)*f(a) > 0:
            a = mid
        else:
            b = mid

        iteration += 1
        error.append((b-a)/2)
        ratio = error[iteration]/error[iteration-1]
        print(f"x_{iteration}: {mid:.6f} ,error: {error[iteration]:.6f} ,e_{iteration}/e_{iteration-1}: {ratio:.6f}")

    return (a+b)/2 ,(b-a)/2

def f(x):
    return x**3

res, err = bisection(f ,-100.565664, 333.5456,1e-10)
print(f"root find:{res}, with error: {err}")

2.1.2 Theoretical Foundations link

For any $n \ge 1$, there exists subintervals $[a_n, b_n]$ of $[a, b]$ such that $b_n - a_n = \frac{1}{2^n} (b-a)$ and $P^* \in (a_n,b_n)$.

$$P_n = \frac{1}{2} (a_n + b_n)$$

$$\therefore |P_n - P^*| \le \frac{1}{2} (b_n - a_n) = \frac{1}{2} \cdot \frac{1}{2^{n-1}} (b-a) = \frac{1}{2^n} (b-a)$$

Hence $\{P_n\}$ converges to $P^*$ with rate $\mathcal{O}\left(\frac{1}{2^n}\right)$.

Note: For a function $f$ with accuracy $\epsilon$ in interval $[a ,b]$, how many steps are required?

$$|P_n - P^*| \le \frac{1}{2^n} (b-a) < \epsilon$$

Where $n$ is the minimum steps required.

2.2 Fixed-Point Iteration link

2.2.1 Definition & Properties link

Explanation:

1. Imagine there is a square with the range $x \in [a ,b]$, $y \in [a, b]$. Then, $g(x)$ is a continuous curve extending from $x = a$ to $x = b$. Because the value of $g(x)$ is always between $[a,b]$, this line must cross the diagonal $y = x$ at least once. That intersection point is the fixed point $X^*$.

2. $|g'(x)|$ is the slope of the function. The restriction $|g'(x)| \le k < 1$ means that the function’s change is relatively “gradual,” even flatter than the diagonal $y = x$. If you cross the diagonal once and want to cross it again, your slope must become steeper than $1$. Since the slope is restricted to below $1$, once the function crosses the diagonal, it can never return.

Existence: Case 1: $g(a) = a$ or $g(b) = b$.

Then $g(x)$ trivially has a fixed point. Case 2: $g(a) > a$ and $g(b) < b$. Let $h(x) := g(x) - x$, which is continuous on $[a,b]$ and $h(a) > 0$, $h(b) < 0$. By the Intermediate Value Theorem, $\exists X^* \in (a ,b)$ such that $h(X^*) = 0$, i.e. $g(X^*) = X^*$. Hence, the fixed point exists.

Uniqueness: Suppose $\exists \ y^*, y^* \neq X^*$ such that $g(y^*) = y^*$. By the Mean Value Theorem, $\exists \ \xi \in (X^*, y^*)$ such that:

$$g'(\xi) = \frac{g(X^*) - g(y^*)}{X^* - y^*}$$

$$\begin{aligned} \therefore |g(X^*) - g(y^*)| &= |g'(\xi)||X^* - y^*| \\ \Rightarrow |X^* - y^*| &= |g'(\xi)||X^* - y^*| \\ \Rightarrow |g'(\xi)| &= 1 \quad (\Rightarrow\Leftarrow) \end{aligned}$$

Contradiction! Hence, $X^* = y^*$, which means the fixed point is unique.

Extension: Why can’t $g'(x)$ be less than $-1$? One might think that even if $g'(x) < -1$, the fixed point is still unique. THIS IS FALSE, because it would never satisfy the self-mapping condition $g(x) \in [a ,b]$. By the Mean Value Theorem, $\exists \ c \in (a ,b)$ such that:

$$\begin{aligned} g(b) - g(a) &= g'(c)(b-a), \quad \text{where } c \in (a, b) \text{ and } g'(x) < -1 \\ &\Rightarrow g(b) - g(a) < -1 \cdot (b-a) \\ &\Rightarrow g(a) - g(b) > b - a \\ &\Rightarrow |g(a) - g(b)| > |a - b| \\ &\therefore g([a,b]) \not\subset [a, b] \end{aligned}$$

2.2.2 Algorithms: Find Fixed Point link

def fixed_point(g, x_0, tolerance = 1e-10, max_iteration = 100):
    x = [x_0]
    for k in range(1, max_iteration+1):
        x.append(g(x[k-1]))
        if abs(x[k] - x[k-1]) < tolerance:
            return x[k]
    return x[-1]

2.2.3 Convergence Theorem & Error Analysis link

For $n > 0$, and $X^*$ as the fixed point:

$$\begin{aligned} |x_n - X^*| &= |g(x_{n-1}) - g(X^*)| \quad \text{(by MVT)} \\ &= |g'(\xi) (x_{n-1} - X^*)|, \quad \xi \in (x_{n-1} ,X^*) \\ &\le k|x_{n-1} - X^*| \le k^2|x_{n-2} - X^*| \le \cdots \le k^n|x_0 - X^*| \end{aligned}$$

Since $k \in (0,1)$, $\lim_{n \rightarrow \infty} k^n = 0 \Rightarrow \lim_{n \rightarrow \infty} |x_n - X^*|=0$. Hence, $\{x_n\}$ converges to $X^*$.

$$\begin{aligned} &\because x_{n+1} - X^* = g(x_n) - g(X^*) = g'(\xi)(x_n - X^*) \quad \text{(MVT)} \\ &\text{Since } g'(x) \text{ is continuous on } (a ,b) \text{, } \lim_{n \rightarrow \infty} g'(\xi) = g'(X^*)\\ &\text{Thus, } \lim_{n \rightarrow \infty} \frac{|x_{n+1} - X^*|}{|x_n - X^*|} = \lim_{n \rightarrow \infty} |g'(\xi)| = |g'(X^*)| \le k \\ &\therefore \alpha = 1 \text{, } \lambda = g'(X^*) \\ &\text{Hence, } \{x_n\} \text{ converges linearly to } X^*. \end{aligned}$$

2.3 Newton’s Method link

If $f \in C^2[a ,b], x_0 \in (a,b)$, the solution $X^* \in (a ,b)$ and $f(X^*) = 0$. Assume $x_0$ is close to $X^*$. By Taylor’s expansion:

$$f(X^*) = f(x_0) + f'(x_0)(X^* - x_0) + \frac{f''(\xi)}{2!} (X^* - x_0)^2, \quad \xi \in (x_0,X^*)$$

$\because |x_0 - X^*|$ is small, $\therefore |x_0 - X^*|^2 \rightarrow 0$ and $\because f \in C^2[a ,b]$, $\therefore |f''(\xi)|$ is bounded.

$$\begin{aligned} 0 &\approx f(x_0) + f'(x_0)(X^* - x_0) \\ \Rightarrow X^* &\approx x_0 - \frac{f(x_0)}{f'(x_0)} \\ \Rightarrow x_1 &= x_0 - \frac{f(x_0)}{f'(x_0)} \\ \Rightarrow x_n &= x_{n-1} - \frac{f(x_{n-1})}{f'(x_{n-1})}, \quad n = 1, 2, \dots \end{aligned}$$

2.3.1 Algorithms link

import math

def d(f, x):
    h = 1e-6
    return (f(x + h) - f(x - h)) / (2 * h)

def f(x):
    return math.log(x) - 1

def newton_method(f, a, tolerance):
    res = [a]
    res.append(res[0] - f(res[0]) / d(f, res[0]))
    iteration = 1

    while (abs(res[iteration] - res[iteration-1]) >= tolerance or abs(f(res[iteration])) >= tolerance):
        current_x = res[iteration]
        next_x = current_x - f(current_x) / d(f, current_x)
        res.append(next_x)
        iteration += 1
        print(f"Iteration {iteration}: x = {res[iteration]:.6f}, f(x) = {f(res[iteration]):.6e}")

    return res[iteration]

final_root = newton_method(f, 4, 1e-6)
print(f"Found root: {final_root:.6f}")

2.3.2 Convergence of Newton’s Method link

Let $g(x) := x - \frac{f(x)}{f'(x)}$. Since $f'$ is continuous and $f'(X^*) \neq 0$, $\exists \bar{\delta} > 0$ such that $f'(x) \neq 0, \forall x \in (X^* - \bar{\delta}, X^* + \bar{\delta})$. Hence, $g(x)$ is continuous on $(X^* - \bar{\delta} ,X^* + \bar{\delta})$ and:

$$g'(x) = 1 - \frac{f'(x)f'(x) - f(x)f''(x)}{[f'(x)]^2} = \frac{f(x)f''(x)}{[f'(x)]^2} \quad \text{and} \quad g'(X^*) = \frac{f(X^*)f''(X^*)}{[f'(X^*)]^2} = 0$$

$$\begin{aligned} \because \ &g'(x) \text{ is continuous around } X^* \text{ and } \exists \ k \in (0, 1) \text{ s.t. } \\ &|g'(x)| \le k < 1, \quad \forall x \in [X^* - \delta_1, X^* + \delta_1], \text{ for some } \delta_1 > 0. \\ &\text{Let } \delta = \min(\delta_1, \bar{\delta}), \text{ then } |g(x) - g(X^*)| = |g'(\xi)||X^* - x| \quad \text{(MVT for $\xi \in (X^*, x)$)} \\ \therefore \ &|g(x) - X^*| = |g(x) - g(X^*)| \le k|x-X^*| < |x-X^*| \end{aligned}$$

Hence, the function $g(x)$ maps $[X^* - \delta_1, X^* + \delta_1]$ onto $[X^* - \delta, X^* + \delta]$. By the fixed point theorem, the sequence $\{x_n\}_{n = 1}^{\infty}$ generated by $x_n = g(x_{n-1}), n \ge 1$ will converge to $X^*$ for any initial value $x_0 \in (X^* - \delta, X^* + \delta)$.

2.3.3 Error Analysis link

3 Interpolation and Polynomial approximation link

3.1 Lagrange Polynomial link

In numerical analysis, when we obtain a set of discrete data points, we often need to construct a polynomial that passes precisely through these points. This process is called Polynomial Interpolation. The Lagrange polynomial provides an intuitive method to construct such a function.

Initial Condition: Suppose we have $n+1$ distinct data points on a function $f(x)$: $x_0, x_1, \dots, x_n$. Our goal is to find a polynomial $P(x)$ of degree at most $n$ that satisfies the interpolation condition for all points: $P(x_i) = f(x_i)$.

3.1.1 Standard Derivation link

First, we construct the Lagrange Basis Polynomials $L_{n,i}(x)$. The key property of this basis function is that it equals $1$ at its corresponding node $x_i$, and $0$ at all other nodes:

$$L_{n,i}(x_j) = \begin{cases} 1 & , \text{ if } i = j \\ 0 & , \text{ if } i \neq j \end{cases}$$

By multiplying the terms $(x-x_k)$ for all other nodes and normalizing the result, we can write the general formula for $L_{n,i}(x)$:

$$L_{n,i}(x) = \frac{(x-x_0)(x-x_1) \cdots (x-x_{i-1})(x-x_{i+1}) \cdots (x-x_n)}{(x_i-x_0)(x_i-x_1) \cdots (x_i-x_{i-1})(x_i-x_{i+1}) \cdots (x_i-x_n)} = \prod_{\substack{k=0 \\ k \neq i}}^{n} \frac{x-x_k}{x_i-x_k}$$

With the basis functions defined, we simply multiply each basis by its corresponding function value $f(x_i)$ and sum them up to obtain the final interpolation polynomial $P(x)$:

$$P(x) = \sum_{i=0}^{n} L_{n,i}(x) \cdot f(x_i) = \sum_{i=0}^{n} \left( \frac{\prod_{j \neq i}^{n} (x-x_j)}{\prod_{j \neq i}^{n} (x_i-x_j)} \right) \cdot f(x_i)$$

3.1.2 Barycentric Form link

The traditional Lagrange polynomial is computationally expensive to evaluate. To resolve this, we introduce an auxiliary Nodal Polynomial $Q(x)$ to unify all the numerator terms:

$$\text{Let } Q(x) = \prod_{j=0}^{n} (x - x_j)$$

To find the relationship between the denominator $\prod_{j \neq i}^{n} (x_i-x_j)$ and $Q(x)$, we apply Logarithmic Differentiation. This technique converts the product into a sum, making differentiation much easier:

$$\begin{aligned} \ln Q(x) &= \sum_{j=0}^{n} \ln(x - x_j) \\ \frac{Q'(x)}{Q(x)} &= \sum_{j=0}^{n} \frac{1}{x - x_j} \quad \text{(Differentiating both sides)} \\ Q'(x) &= \sum_{j=0}^{n} \frac{Q(x)}{x - x_j} \end{aligned}$$

Next, we evaluate $Q'(x)$ precisely at $x = x_i$. Since $Q(x_i) = 0$, all terms in the summation evaluate to $0$ except for the term where the denominator is $(x_i - x_i)$, which is preserved by evaluating the limit (L’Hôpital’s rule):

$$Q'(x_i) = \lim_{x \to x_i} \sum_{j=0}^{n} \frac{Q(x)}{x - x_j} = \prod_{\substack{j=0 \\ j \neq i}}^{n} (x_i - x_j)$$

This result is exactly the denominator in our original $L_{n,i}(x)$ formula! Therefore, we can redefine a much more computationally efficient basis function $\phi_{n,i}(x)$:

$$\phi_{n,i}(x) = \frac{\frac{Q(x)}{x-x_i}}{Q'(x_i)} = \frac{Q(x)}{(x-x_i)Q'(x_i)}$$

Finally, the Lagrange polynomial can be elegantly rewritten into the First Form of Barycentric Interpolation:

$$P(x) = \sum_{i=0}^{n} \frac{Q(x)}{(x-x_i)Q'(x_i)} \cdot f(x_i)$$

This form drastically reduces computational complexity. If a new data point is added, we only need to update the weights rather than recalculating all polynomials from scratch.

3.2 Error Analysis of Lagrange Interpolation link

Case 1: $x = x_k$ If $x = x_k$ for any $k = 0, \dots, n$, then $f(x_k) = P(x_k)$, and $\xi(x)$ can be arbitrarily chosen in $(a,b)$.

Case 2: $x \neq x_k$ If $x \neq x_k$ for all $k = 0, \dots, n$, we define a function $g$ for $t \in [a,b]$:

$$g(t) = f(t) - P(t) - [f(x) - P(x)] \prod_{i=0}^{n} \frac{t-x_i}{x-x_i}$$

Let $t = x$, we have:

$$g(x) = f(x) - P(x) - [f(x) - P(x)] \prod_{i=0}^{n} \frac{x-x_i}{x-x_i} = f(x) - P(x) - [f(x) - P(x)] = 0$$

Thus $g \in C^{n+1}[a,b]$, and $g(t) = 0$ for all $t \in \{x, x_0, \dots, x_n\}$. By Rolle’s Theorem, there exists $\xi \in (a,b)$ such that $g^{(n+1)}(\xi) = 0$, i.e.,

$$\begin{aligned} 0 = g^{(n+1)}(\xi) &= f^{(n+1)}(\xi) - P^{(n+1)}(\xi) - [f(x) - P(x)] \frac{d^{n+1}}{dt^{n+1}} \left[ \prod_{i=0}^{n} \frac{t-x_i}{x-x_i} \right]_{t=\xi} \\ &= f^{(n+1)}(\xi) - 0 - [f(x) - P(x)] \frac{d^{n+1}}{dt^{n+1}} \left[ \frac{t^{n+1}}{\prod_{i=0}^{n} (x-x_i)} \right]_{t=\xi} \\ &= f^{(n+1)}(\xi) - [f(x) - P(x)] \frac{(n+1)!}{\prod_{i=0}^{n} (x-x_i)} \end{aligned}$$

$$\Rightarrow f(x) = P(x) + \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^{n} (x-x_i)$$

3.3 Newton’s Divided Difference link

Suppose $P_n(x)$ is the $n$-th Lagrange interpolating polynomial with $f(x)$ at distinct numbers $x_0, \dots, x_n$. Although the polynomial is unique, we can alternate its representation. Let

$$P_n(x) = a_0 + a_1(x-x_0) + \dots + a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1})$$

Then, $a_0 = P(x_0) = f(x_0)$, and $a_1 = \frac{f(x_1) - f(x_0)}{x_1 - x_0}$.

To evaluate $a_k$ for $k = 0, \dots, n$, we define the divided difference:

1. $f[x_i] = f(x_i)$

2. $f[x_i, x_{i+1}] = \frac{f[x_{i+1}] - f[x_i]}{x_{i+1} - x_i}$

3. $f[x_1, \dots, x_n] = \frac{f[x_2, \dots, x_n] - f[x_1, \dots, x_{n-1}]}{x_n - x_1}$

Then we can write $f[x_0] = a_0$, $f[x_0, x_1] = a_1$. As might be expected from evaluating $a_0$ and $a_1$, $a_k = f[x_0, \dots, x_k]$. Thus, $P_n(x)$ can be rewritten in a form called Newton’s Divided Differences:

$$P_n(x) = f[x_0] + \sum_{k=1}^{n} f[x_0, \dots, x_k] (x-x_0)(x-x_1)\cdots(x-x_{k-1})$$

Let $g(x) = f(x) - P_n(x)$. Since $f(x_i) = P_n(x_i)$ for $i = 0, \dots, n$, $g$ has $n+1$ distinct zeros in $[a,b]$. By Rolle’s Theorem, there exists $\xi \in (a,b)$ such that $g^{(n)}(\xi) = 0$.

$$\Rightarrow f^{(n)}(\xi) - P_n^{(n)}(\xi) = 0 \Rightarrow f^{(n)}(\xi) = P_n^{(n)}(\xi)$$

Since $P_n(x) = f[x_0, \dots, x_n]x^n + \dots$, it follows that $P_n^{(n)}(x) = n! f[x_0, \dots, x_n]$. Hence, $f^{(n)}(\xi) = n! f[x_0, \dots, x_n]$.

If the nodes are arranged uniformly, Newton’s Divided Difference formula can be expressed in a simplified form. Assume $x_i = x_0 + ih$ and let $x = x_0 + sh$, the difference becomes $x - x_i = (s-i)h$. So,

$$\begin{aligned} P_n(x) &= P_n(x_0 + sh) \\ &= f[x_0] + sh f[x_0, x_1] + sh(s-1)h f[x_0, x_1, x_2] + \dots + s(s-1)\cdots(s-k+1)h^k f[x_0, \dots, x_k] + \dots \\ &= f(x_0) + \sum_{k=1}^{n} \binom{s}{k} k! h^k f[x_0, \dots, x_k] \\ &= f(x_0) + \sum_{k=1}^{n} \binom{s}{k} \Delta^k f(x_0) \end{aligned}$$

3.4 Cubic Spline Interpolation link

To construct the cubic spline interpolant for a function $f$:

$$S_j(x) = a_j + b_j(x-x_j) + c_j(x-x_j)^2 + d_j(x-x_j)^3$$

Since $a_j = S_j(x_j)$, and by definition condition (3):

$$a_{j+1} = S_{j+1}(x_{j+1}) = S_j(x_{j+1}) = a_j + b_j(x_{j+1}-x_j) + c_j(x_{j+1}-x_j)^2 + d_j(x_{j+1}-x_j)^3$$

Let $h_j = x_{j+1} - x_j$ for $j = 0, \dots, n-2$, then:

$$a_{j+1} = a_j + b_j h_j + c_j h_j^2 + d_j h_j^3 \quad \text{for } j = 0, \dots, n-1 \label{eq:spline1}$$

Similarly, from $b_j = S'_j(x_j)$, we have:

$$b_{j+1} = b_j + 2c_j h_j + 3d_j h_j^2 \quad \text{for } j = 0, \dots, n-1 \label{eq:spline2}$$

Similarly, from $2c_j = S''_j(x_j)$, we have $2c_{j+1} = 2c_j + 6d_j h_j$, which yields:

$$d_j = \frac{c_{j+1} - c_j}{3h_j} \label{eq:spline4}$$

Substitute [eq:spline4] into [eq:spline1] and [eq:spline2]:

$$a_{j+1} = a_j + b_j h_j + \frac{h_j^2}{3}(2c_j + c_{j+1})$$

$$b_{j+1} = b_j + h_j(c_j + c_{j+1}) \label{eq:spline5}$$

Isolating $b_j$ and $b_{j-1}$:

$$\begin{aligned} b_j &= \frac{1}{h_j}(a_{j+1} - a_j) - \frac{h_j}{3}(2c_j + c_{j+1}) \\ b_{j-1} &= \frac{1}{h_{j-1}}(a_j - a_{j-1}) - \frac{h_{j-1}}{3}(2c_{j-1} + c_j) \end{aligned}$$

Substitute $b_j$ and $b_{j-1}$ (as $b_{j+1}$ and $b_j$) into [eq:spline5], we get:

$$h_{j-1} c_{j-1} + 2(h_{j-1} + h_j)c_j + h_j c_{j+1} = \frac{3}{h_j}(a_{j+1} - a_j) - \frac{3}{h_{j-1}}(a_j - a_{j-1})$$

Where $\{c_j\}_{j=0}^{n}$ are unknowns, $\{h_j\}_{j=0}^{n-1}$ and $\{a_j\}_{j=0}^{n}$ are given, and $\{b_j\}_{j=0}^{n-1}$ and $\{d_j\}_{j=0}^{n-1}$ depend on $c_j$. Thus, once $\{c_j\}_{j=0}^{n}$ are determined, we can construct $\{S_j(x)\}_{j=0}^{n-1}$.

To find $\{c_j\}_{j=0}^{n}$ using the definition $S''(x_0) = S''(x_n) = 0$, i.e., $c_0 = c_n = 0$. From the equation above, we can construct a linear system $Ax = b$:

$$A = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ h_0 & 2(h_0+h_1) & h_1 & \cdots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1} \\ 0 & \cdots & 0 & 0 & 1 \end{bmatrix} \in \mathbb{R}^{(n+1) \times (n+1)}$$

$$x = \begin{bmatrix} c_0 \\ \vdots \\ c_n \end{bmatrix}, \quad b = \begin{bmatrix} 0 \\ \frac{3}{h_1}(a_2-a_1) - \frac{3}{h_0}(a_1-a_0) \\ \vdots \\ \frac{3}{h_{n-1}}(a_n-a_{n-1}) - \frac{3}{h_{n-2}}(a_{n-1}-a_{n-2}) \\ 0 \end{bmatrix}$$

Since $A$ is strictly diagonally dominant ($|a_{ii}| > \sum_{j \neq i}^{n+1} |a_{ij}|$), the system has a unique solution.

$\because S''(x_0) = S''(x_n) = 0$ $\therefore c_n = S''(x_n)/2 = 0$ $S''(x_0) = 2c_0 + 6d_0(x_0 - x_0) = 0 \Rightarrow c_0 = 0$ By deriving $c_0 = c_n$, we can construct a linear system: $Ax = b$, where $A, x, b$ is the same as above. Since $A$ is strictly diagonally dominant, the system has a unique solution.

4 Numerical Differentiation and Integration link

4.1 Numerical Differentiation link

To evaluate $f'(x)$ at $x = x_0$, we can simply compute:

$$f'(x_0) \approx \frac{f(x_0+h) - f(x_0)}{h}$$

for a small value $h$. But it is not always success due to the round error in computer. Hence, we construct the first Lagrange polynomial $P_{0,1}(x)$ for $f$ determined by $x_0, x_1$.

$$f(x) = P_{0,1}(x) + \frac{(x-x_0)(x-x_1)}{2!} f''(\xi(x))$$

$$= f(x_0) \frac{x-x_1}{x_0-x_1} + f(x_1) \frac{x-x_0}{x_1-x_0} + \frac{(x-x_0)(x-x_1)}{2!} f''(\xi(x))$$

Since $x_1 = x_0 + h$, for some $\xi(x) \in (x_0, x_1)$:

$$\begin{aligned} f'(x) &= \frac{f(x_0+h) - f(x_0)}{h} + D_x \left[ \frac{(x-x_0)(x-x_1)}{2!} f''(\xi(x)) \right] \\ &= \frac{f(x_0+h) - f(x_0)}{h} + \frac{2(x-x_0)-h}{2} f''(\xi(x)) + \frac{(x-x_0)(x-x_1)}{2!} D_x[f''(\xi(x))] \end{aligned}$$

Take $x = x_0$:

$$f'(x_0) = \frac{f(x_0+h) - f(x_0)}{h} - \frac{h}{2} f''(\xi(x_0))$$

Suppose $M = \max(|f''(x)|)$ for $x \in (x_0, x_1)$, then we have Error $E \triangleq \left|\frac{h}{2} f''(\xi(x))\right| \le \frac{h}{2} M$. This formula is known as forward-difference formula if $h > 0$, backward-difference formula if $h < 0$.

To obtain general derivates approximation formula, suppose $\{x_i\}_{i=0}^{n}$, $x_i \neq x_j \ \forall i \neq j$ in some interval $I$ and $f \in C^{n+1}(I)$. By Lagrange polynomial interpolation:

$$f(x) = \sum_{k=0}^{n} f(x_k) L_k(x) + \frac{\prod_{j=0}^{n} (x-x_j)}{(n+1)!} f^{(n+1)}(\xi(x))$$

for some $\xi(x) \in I$. $L_k(x)$ is the $k$-th Lagrange polynomial. Differentiating this expression:

$$f'(x) = \sum_{k=0}^{n} f(x_k) L'_k(x) + D_x \left[ \frac{\prod_{j=0}^{n} (x-x_j)}{(n+1)!} \right] f^{(n+1)}(\xi(x)) + \frac{\prod_{j=0}^{n} (x-x_j)}{(n+1)!} D_x [f^{(n+1)}(\xi(x))]$$

Let $x = x_i$, $i = 0 \dots n$:

$$f'(x_i) = \sum_{k=0}^{n} f(x_k) L'_k(x_i) + \frac{f^{(n+1)}(\xi(x_i))}{(n+1)!} \prod_{\substack{j=0 \\ j \neq i}}^{n} (x_i - x_j)$$

This formula is called $(n+1)$-point formula. Usually, we use $f'(x_i) \approx \sum_{k=0}^{n} f(x_k) L'_k(x_i)$.

4.1.1 Three-Points Formulas link

Suppose $x_1 = x_0 + h$, $x_2 = x_0 + 2h$, $h \neq 0$:

$$\begin{aligned} L_0(x) &= \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}, \quad L'_0(x) = \frac{1}{2h^2}(2x-x_1-x_2) \\ L_1(x) &= \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}, \quad L'_1(x) = -\frac{1}{h^2}(2x-x_0-x_2) \\ L_2(x) &= \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}, \quad L'_2(x) = \frac{1}{2h^2}(2x-x_0-x_1) \end{aligned}$$

$$\begin{aligned} f'(x_j) &= f(x_0)\left[ \frac{1}{2h^2}(2x_j-x_1-x_2) \right] + f(x_1)\left[ -\frac{1}{h^2}(2x_j-x_0-x_2) \right] \\ &\quad + f(x_2)\left[ \frac{1}{2h^2}(2x_j-x_0-x_1) \right] + \frac{f^{(3)}(\xi_j)}{3!} \prod_{\substack{k=0 \\ k \neq j}}^{2} (x_j - x_k) \end{aligned}$$

Let $x_j = x_0$:

$$f'(x_0) = \frac{1}{h} \left[ -\frac{3}{2}f(x_0) + 2f(x_1) - \frac{1}{2}f(x_2) \right] + \frac{h^2}{3} f^{(3)}(\xi_0)$$

Let $x_j = x_1$:

$$f'(x_1) = \frac{1}{h} \left[ -\frac{1}{2}f(x_0) + \frac{1}{2}f(x_2) \right] - \frac{h^2}{6} f^{(3)}(\xi_1)$$

Let $x_j = x_2$:

$$f'(x_2) = \frac{1}{h} \left[ \frac{1}{2}f(x_0) - 2f(x_1) + \frac{3}{2}f(x_2) \right] + \frac{h^2}{3} f^{(3)}(\xi_2)$$

Substitute $x_1 = x_0 + h$, $x_2 = x_0 + 2h$ can get another expression.

4.1.2 Five-Points Formulas link

$$f'(x_0) = \frac{1}{12h} [-25f(x_0) + 48f(x_0+h) - 36f(x_0+2h) + 16f(x_0+3h) - 3f(x_0+4h)] + \frac{h^4}{5} f^{(5)}(\xi), \quad \xi \in (x_0, x_0+4h)$$

4.1.3 Second Derivative Midpoint Formula link

Consider the approximation for $f''(x_0)$ and $f \in C^4[x_0-h, x_0+h]$. The most direct method is to use Taylor expansions on $f(x_0+h)$ and $f(x_0-h)$:

$$\begin{aligned} f(x_0+h) &= f(x_0) + f'(x_0)h + \frac{h^2}{2!}f''(x_0) + \frac{h^3}{3!}f'''(x_0) + \frac{h^4}{4!}f^{(4)}(\xi_1) \label{eq:taylor1} \\ f(x_0-h) &= f(x_0) - f'(x_0)h + \frac{h^2}{2!}f''(x_0) - \frac{h^3}{3!}f'''(x_0) + \frac{h^4}{4!}f^{(4)}(\xi_2) \label{eq:taylor2} \end{aligned}$$

By adding [eq:taylor1] and [eq:taylor2], we can cancel out the terms with odd powers of $h$, then we get:

$$f(x_0+h) + f(x_0-h) = 2f(x_0) + h^2 f''(x_0) + \frac{h^4}{24}[f^{(4)}(\xi_1) + f^{(4)}(\xi_2)]$$

where $\xi_1 \in (x_0, x_0+h)$, $\xi_2 \in (x_0-h, x_0)$. Since $f \in C^4[x_0-h, x_0+h]$ and $f^{(4)}$ is bounded on the interval. By IVT, $\exists \xi \in [x_0-h, x_0+h]$ s.t. $\frac{1}{2}(f^{(4)}(\xi_1) + f^{(4)}(\xi_2)) = f^{(4)}(\xi)$. Substitute it back, we have:

$$f(x_0+h) + f(x_0-h) = 2f(x_0) + h^2 f''(x_0) + \frac{h^4}{12} f^{(4)}(\xi)$$

$$\Rightarrow f''(x_0) = \frac{1}{h^2}[f(x_0-h) - 2f(x_0) + f(x_0+h)] - \frac{h^2}{12}f^{(4)}(\xi)$$

4.1.4 Round-Off Error Instability link

Particularly, it is necessary to pay attention to round-off error. Suppose our computation actually use $\tilde{f}(x_0+h)$ and $\tilde{f}(x_0-h)$ with round-off error $e(x_0+h)$ and $e(x_0-h)$ respectively. i.e.

$$\begin{cases} f(x_0+h) = \tilde{f}(x_0+h) + e(x_0+h) \\ f(x_0-h) = \tilde{f}(x_0-h) + e(x_0-h) \end{cases}$$

If we use these terms to approximate $f'(x_0)$, the total error

$$E = f'(x_0) - \frac{\tilde{f}(x_0+h) - \tilde{f}(x_0-h)}{2h} = \frac{1}{2h}[e(x_0+h) - e(x_0-h)] - \frac{h^2}{6}f^{(3)}(\xi_1)$$

Suppose $|e(x_0 \pm h)| \le \varepsilon$ and $|f^{(3)}(\xi_1)| \le M$, then

$$|E| \le \frac{\varepsilon}{h} + \frac{h^2}{6}M$$

If $\varepsilon$ is given, $E$ is a function only regarding to $h$. We can find the best $h$ to minimize the error:

$$E'(h) = -\frac{\varepsilon}{h^2} + \frac{h}{3}M = 0 \implies h = \sqrt[3]{\frac{3\varepsilon}{M}}$$

4.2 Numerical Integration link

4.2.1 The Trapezoidal Rule link

To derive the Trapezoidal Rule for approximating $\int_a^b f(x) dx$, let $a = x_0$, $b = x_1$, $h = b - a$, and we use Lagrange polynomial.

$$P_1(x) = \frac{(x-x_1)}{(x_0-x_1)}f(x_0) + \frac{(x-x_0)}{(x_1-x_0)}f(x_1)$$

then,

$$\int_a^b f(x) dx = \int_{x_0}^{x_1} P_1(x) + E(x) dx$$

$$\begin{aligned} \int_{x_0}^{x_1} P_1(x) dx &= \int_{x_0}^{x_1} \left[ f(x_0)\frac{x-x_1}{x_0-x_1} + f(x_1)\frac{x-x_0}{x_1-x_0} \right] dx \\ &= \left[ f(x_0) \frac{(x-x_1)^2}{2(x_0-x_1)} + f(x_1) \frac{(x-x_0)^2}{2(x_1-x_0)} \right]_{x_0}^{x_1} \\ &= \frac{(x_1-x_0)}{2}[f(x_0) + f(x_1)] \end{aligned}$$

$$\int_{x_0}^{x_1} E(x) dx = \frac{1}{2} \int_{x_0}^{x_1} f''(\xi(x))(x-x_0)(x-x_1) dx$$

$\because f \in C^2[a,b]$. By Extreme Value Theorem, $m < f''(x) < M \ \forall x \in [x_0, x_1]$. Let $g(x) = (x-x_0)(x-x_1)$, $E(x) = \frac{1}{2}f''(x)g(x)$ i.e. $m \cdot g(x) \le E(x) \le M \cdot g(x) \Rightarrow \int_{x_0}^{x_1} m \, g(x) dx \le \int_{x_0}^{x_1} E(x) dx \le \int_{x_0}^{x_1} M \, g(x) dx$. Due to $m$ and $M$ are constants, we have

$$m \int_{x_0}^{x_1} g(x) dx \le \int_{x_0}^{x_1} E(x) dx \le M \int_{x_0}^{x_1} g(x) dx$$

$$\Rightarrow m \le \frac{\int_{x_0}^{x_1} E(x) dx}{\int_{x_0}^{x_1} g(x) dx} \le M$$

By IVT, $\exists \xi \in (x_0, x_1)$ s.t.

$$\int_{x_0}^{x_1} E(x) dx = \frac{1}{2} f''(\xi) \int_{x_0}^{x_1} (x-x_0)(x-x_1) dx = -\frac{h^3}{12} f''(\xi)$$

Thus, substitute $P_1(x)$ and $E(x)$:

$$\int_{x_0}^{x_1} f(x) dx = \frac{h}{2} [f(x_0) + f(x_1)] - \frac{h^3}{12} f''(\xi)$$

4.2.2 The Simpson’s Rule link

Simpson’s Rule uses one more point for Lagrange polynomial. $x_0 = a$, $x_1 = a+h$, $x_2 = b$, $h = \frac{b-a}{2}$.

$$\begin{aligned} \int_a^b f(x) dx &= \int_{x_0}^{x_2} f(x_0) \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} dx + \int_{x_0}^{x_2} f(x_1) \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} dx \\ &\quad + \int_{x_0}^{x_2} f(x_2) \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} dx + \int_{x_0}^{x_2} f^{(3)}(\xi(x)) (x-x_0)(x-x_1)(x-x_2) \cdot \frac{1}{3!} dx \\ &= \frac{h}{3} f(x_0) + \frac{4h}{3} f(x_1) + \frac{h}{3} f(x_2) + E(x) \end{aligned}$$

$$E(x) = \frac{1}{6} \int_{x_0}^{x_2} f^{(3)}(\xi(x)) (x-x_0)(x-x_1)(x-x_2) dx$$

Due to $E(x)$ now is an odd function, $E(x) = 0$. Thus, the error is not determined on a polynomial with $3$ degree. Alternatively, we use Hermite interpolation by adding a restriction $P_3'(x_1) = f'(x_1)$.

$$\begin{aligned} E(x) &= \frac{1}{4!} \int_{x_0}^{x_2} f^{(4)}(\xi(x)) (x-x_0)(x-x_1)^2(x-x_2) dx \\ &= \frac{1}{24} f^{(4)}(\xi) \int_{x_0}^{x_2} (x-x_0)(x-x_1)^2(x-x_2) dx \quad \text{(WMVT)}, \ \xi \in (a,b) \\ &= \frac{1}{24} f^{(4)}(\xi) \int_{-h}^{h} (t+h)(t^2)(t-h) dt \quad (t = x-x_1) \\ &= \frac{1}{24} f^{(4)}(\xi) \int_{-h}^{h} t^4 - h^2 t^2 dt = -\frac{h^5}{90} f^{(4)}(\xi) \end{aligned}$$

4.3 Adaptive Quadrature Methods link

Adaptive Quadrature is to logically optimize computational resources. Standard composite formulas use a uniform step size, which oversamples regions of low functional variation and undersamples regions of high variation.

Let $S(a,b)$ be the Simpson’s Rule for $\int_a^b f(x) dx$. $S(a,b) \triangleq \frac{h}{3}[f(a) + 4f(a+h) + f(b)]$, where $h = \frac{b-a}{2}$, then

$$\begin{cases} S(a, \frac{a+b}{2}) = \frac{h}{6}[f(a) + 4f(a+\frac{h}{2}) + f(a+h)] \\ S(\frac{a+b}{2}, b) = \frac{h}{6}[f(a+h) + 4f(a+\frac{3}{2}h) + f(a+2h)] \end{cases}$$

$$\begin{aligned} \int_a^b f(x) dx &= S(a, \frac{a+b}{2}) + S(\frac{a+b}{2}, b) - \left(\frac{1}{2}\right)^4 \frac{h^5}{90} f^{(4)}(\xi), \quad \xi \in (a,b) \\ &= S(a, b) - \frac{h^5}{90} f^{(4)}(\xi_1), \quad \xi_1 \in (a,b) \end{aligned}$$

Suppose $\xi_1 \approx \xi$, $f^{(4)}(\xi_1) \approx f^{(4)}(\xi)$:

$$S(a, \frac{a+b}{2}) + S(\frac{a+b}{2}, b) - \frac{1}{16} \frac{h^5}{90} f^{(4)}(\xi) \approx S(a, b) - \frac{h^5}{90} f^{(4)}(\xi_1)$$

$$\Rightarrow \underbrace{\frac{h^5}{90} f^{(4)}(\xi)}_{\text{Error}} \approx \frac{16}{15} \left[ S(a, b) - S(a, \frac{a+b}{2}) - S(\frac{a+b}{2}, b) \right]$$

4.4 Gaussian Quadrature link

Gaussian quadrature chooses the points for evaluation in an optimal way. The $x_1 \dots x_n$ in the interval, $c_1 \dots c_n$ are chosen to minimize the expected error obtained in the approximation.

$$\int_a^b f(x) dx \approx \sum_{i=1}^n c_i f(x_i), \quad x_i \in [a,b] \ \forall i=1\sim n$$

Consider $n=2$, $[a,b] = [-1,1]$:

$$\int_{-1}^1 f(x) dx \approx c_1 f(x_1) + c_2 f(x_2)$$

Due to we have four parameters, we construct a 3-order polynomial to determine these parameter. Then,

$$f(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3 \approx c_1 f(x_1) + c_2 f(x_2)$$

$$\int_{-1}^1 f(x) dx = a_0 \int_{-1}^1 1 dx + a_1 \int_{-1}^1 x dx + a_2 \int_{-1}^1 x^2 dx + a_3 \int_{-1}^1 x^3 dx$$

If

$$\begin{aligned} f(x)=1 &\Rightarrow c_1 + c_2 = 2 \\ f(x)=x &\Rightarrow c_1 x_1 + c_2 x_2 = 0 \\ f(x)=x^2 &\Rightarrow c_1 x_1^2 + c_2 x_2^2 = \frac{2}{3} \\ f(x)=x^3 &\Rightarrow c_1 x_1^3 + c_2 x_2^3 = 0 \end{aligned}$$

$$\Rightarrow \begin{cases} c_1 = 1 \\ c_2 = 1 \\ x_1 = -\frac{\sqrt{3}}{3} \\ x_2 = \frac{\sqrt{3}}{3} \end{cases}$$

which gives the formula $\int_{-1}^1 f(x) dx = f\left(-\frac{\sqrt{3}}{3}\right) + f\left(\frac{\sqrt{3}}{3}\right)$.

4.4.1 Legendre Polynomials link

This polynomial can get the parameter needed in Gaussian quadrature more easily. Legendre Polynomials is a collection $\{P_0(x), P_1(x), \dots, P_n(x)\}$ with the following properties:

1. For each $n$, $P_n(x)$ is a monic polynomial of degree $n$ ($a_n = 1$).

2. $\int_{-1}^1 P(x) P_n(x) dx = 0$, whenever $P(x)$ is a polynomial $\deg(P(x)) < n$ which means $P_n(x)$ and $P(x)$ are orthogonal to each other.

Case 1: $\deg(f(x)) < n$

$$f(x) = \sum_{i=1}^n L_i(x) f(x_i)$$

and

$$\begin{aligned} \int_{-1}^1 f(x) dx &= \int_{-1}^1 \sum_{i=1}^n L_i(x) f(x_i) dx \\ &= \sum_{i=1}^n \int_{-1}^1 L_i(x) f(x_i) dx \\ &= \sum_{i=1}^n \left( \int_{-1}^1 L_i(x) dx \right) f(x_i) \\ &= \sum_{i=1}^n c_i f(x_i) \end{aligned}$$

Case 2: $n \le \deg(f(x)) < 2n$ Let $P(x) = Q(x)P_n(x) + R(x)$, where $Q(x)$, $R(x)$ are polynomial with degree less than $n$, $P_n(x)$ is the Legendre polynomial with order $n$, $x_i$ are the roots of $P_n(x)$. Thus for $i=1\dots n$:

$$f(x_i) = Q(x_i)P_n(x_i) + R(x_i) = R(x_i)$$

$$\int_{-1}^1 f(x) dx = \int_{-1}^1 R(x) dx = \sum_{i=1}^n c_i R(x_i)$$

Consider evaluating $\int_a^b f(x) dx$ over an arbitrary $[a,b]$. Let $t = \frac{2x-a-b}{b-a}$, $x = \frac{1}{2}[(b-a)t + a + b]$, $dt = \frac{2}{b-a} dx$:

$$\int_a^b f(x) dx = \int_{-1}^1 f\left(\frac{(b-a)t + (b+a)}{2}\right) \cdot \frac{(b-a)}{2} dt$$

4.5 Improper Integrals link

4.5.1 Left Endpoint Singularity link

Consider $\int_a^b \frac{g(x)}{(x-a)^p} dx$, which converges iff $0 < p < 1$, where $g(x) \in C^4[a,b]$ and $0 < p < 1$. In this case, we can construct the forth Taylor Polynomial $P_4(x)$ at $x=a$ of $g(x)$.

$$P_4(x) = g(a) + g'(a)(x-a) + \frac{g''(a)}{2!}(x-a)^2 + \frac{g'''(a)}{3!}(x-a)^3 + \frac{g^{(4)}(a)}{4!}(x-a)^4$$

$\because g(x) = P_4(x) + R_4(x) \Rightarrow R_4(x) = g(x) - P_4(x)$

$$\therefore f(x) = \frac{g(x)}{(x-a)^p} = \frac{P_4(x)}{(x-a)^p} + \frac{R_4(x)}{(x-a)^p} = \frac{P_4(x)}{(x-a)^p} + \frac{g(x) - P_4(x)}{(x-a)^p}$$

$$\int_a^b \frac{P_4(x)}{(x-a)^p} dx = \sum_{k=0}^4 \frac{g^{(k)}(a)}{k!(k+1-p)} (b-a)^{k+1-p}$$

To determine $\int_a^b \frac{g(x)-P_4(x)}{(x-a)^p} dx$, Define

$$G(x) = \begin{cases} \frac{g(x)-P_4(x)}{(x-a)^p} & \text{if } a < x \le b \\ 0 & \text{if } x = a \end{cases}$$

This gives $G(x)$ continuous on $[a,b]$.

5 Solve the Linear System of Equations link

There are two methods to solve the linear system of equations:

• Direct Methods: These methods find the exact solution in a finite number of steps. Examples include Gaussian elimination, LU decomposition, and Cholesky decomposition.

• Iterative Methods: These methods start with an initial guess and iteratively refine it to approach the solution. Examples include Jacobi method, Gauss-Seidel method, and Conjugate Gradient method.

5.1 Matrix Norms and Condition Number link

To effectively evaluate the stability of a linear system and mathematically bound the errors generated by these solvers, we must explicitly define vector and matrix norms.

5.1.1 Vector Norms link

For a vector $x = [x_1, x_2, \dots, x_n]^T \in \mathbb{R}^n$:

• 1-norm: $\|x\|_1 \triangleq \sum_{i=1}^n |x_i|$.

• p-norm ($p>0$): $\|x\|_p \triangleq \left(\sum_{i=1}^n |x_i|^p\right)^{\frac{1}{p}}$.

• Infinity-norm: $\|x\|_\infty \triangleq \max_{1 \le i \le n} |x_i|$.

5.1.2 Matrix Norms link

For a matrix $A \in \mathbb{R}^{m \times n}$ with entries $a_{ij}$:

• Frobenius Norm (F-norm): $\|A\|_F \triangleq \left(\sum_{i=1}^m \sum_{j=1}^n a_{ij}^2\right)^{\frac{1}{2}}$.

• Matrix p-norm: $\|A\|_p \triangleq \sup_{x \in \mathbb{R}^n, x \neq 0} \frac{\|Ax\|_p}{\|x\|_p}$.

• Column sum norm (1-norm): $\|A\|_1 = \max_{1 \le j \le n} \sum_{i=1}^m |a_{ij}|$.

• Row sum norm ($\infty$-norm): $\|A\|_\infty = \max_{1 \le i \le m} \sum_{j=1}^n |a_{ij}|$.

5.1.3 Condition Number and Error Bounds link

Consider a nonsingular matrix $A$. If $x^*$ is the exact analytical solution to $Ax^* = b$ and $x^k$ is a numerical approximation, we strictly define the residual vector as $r = b - Ax^k$.

By definition, $r = b - Ax^k = Ax^* - Ax^k = A(x^* - x^k)$.

$$\Rightarrow x^* - x^k = A^{-1}r \Rightarrow \|x^* - x^k\| \le \|A^{-1}\| \|r\|$$

Additionally, since $b = Ax^*$, we have $\|b\| \le \|A\| \|x^*\|$, which implies $\|x^*\| \ge \frac{\|b\|}{\|A\|}$. Combining these inequalities gives:

$$\frac{\|x^* - x^k\|}{\|x^*\|} \le \frac{\|A^{-1}\| \|r\|}{\|b\| / \|A\|} = \|A\| \|A^{-1}\| \frac{\|r\|}{\|b\|}$$

If $\kappa(A) \approx 1$, the matrix is well-conditioned. If $\kappa(A) \gg 1$, the matrix is considered ill-conditioned, and numerical solutions will be highly sensitive to perturbations.

5.2 Direct Methods link

5.2.1 Gaussian Elimination link

The structural process of Gaussian elimination can be mathematically modeled as a sequential application of elementary linear transformations. Consider a 2-dimensional vector $x = \begin{bmatrix} x_1 & x_2 \end{bmatrix}^T$ where the pivot $x_1 \neq 0$. There exists an elementary transformation matrix $M$ such that the second entry is eliminated:

$$Mx = \begin{bmatrix} x_1 \\ 0 \end{bmatrix}$$

By defining the multiplier $\alpha = -x_2/x_1$, the transformation matrix $M$ takes the form:

$$M = \begin{bmatrix} 1 & 0 \\ \alpha & 1 \end{bmatrix}$$

More generally, in $\mathbb{R}^n$, let $x = \begin{bmatrix} x_1 & \cdots & x_n \end{bmatrix}^T$ with the $k$-th pivot $x_k \neq 0$. We seek an elimination matrix $M_k$ that annihilates all entries strictly below the pivot:

$$M_k x = \begin{bmatrix} x_1 \\ \vdots \\ x_k \\ 0 \\ \vdots \\ 0 \end{bmatrix}$$

This transformation is defined algebraically as a rank-1 update to the identity matrix:

$$M_k = I - \alpha e_k^T$$

where $e_k$ is the $k$-th standard basis vector (indicating the $k$-th position), and the vector $\alpha$ contains the essential multipliers:

$$\alpha = \begin{bmatrix} 0 \\ \vdots \\ 0 \\ \alpha_{k+1} \\ \vdots \\ \alpha_n \end{bmatrix}, \quad \text{with } \alpha_i = \frac{x_i}{x_k} \quad \text{for } k+1 \le i \le n$$

Remark: The negative sign for elimination is inherently absorbed by the subtraction operator in the definition $I - \alpha e_k^T$.

By systematically applying these elementary lower triangular matrices (Frobenius matrices) to a dense matrix $A$, we execute the forward elimination phase. Consider a generic $4 \times 4$ dense matrix $A$. The geometric progression of zeroing out sub-diagonal entries is visualized as follows:

$$A = \begin{bmatrix} \times & \times & \times & \times \\ \times & \times & \times & \times \\ \times & \times & \times & \times \\ \times & \times & \times & \times \end{bmatrix} \quad \xrightarrow{M_1} \quad M_1 A = \begin{bmatrix} \times & \times & \times & \times \\ 0 & \otimes & \otimes & \otimes \\ 0 & \otimes & \otimes & \otimes \\ 0 & \otimes & \otimes & \otimes \end{bmatrix}$$$$\xrightarrow{M_2} \quad M_2 M_1 A = \begin{bmatrix} \times & \times & \times & \times \\ 0 & \otimes & \otimes & \otimes \\ 0 & 0 & \star & \star \\ 0 & 0 & \star & \star \end{bmatrix} \quad \xrightarrow{M_3} \quad M_3 M_2 M_1 A = \begin{bmatrix} \times & \times & \times & \times \\ 0 & \otimes & \otimes & \otimes \\ 0 & 0 & \star & \star \\ 0 & 0 & 0 & \uparrow \end{bmatrix}$$

The conclusive matrix product, $M_3 M_2 M_1 A$, successfully reduces the original system into an upper triangular matrix. This state mathematically guarantees that the exact solution can be swiftly extracted via backward substitution.

5.2.2 Algorithm link

Given $A \in \mathbb{R}^{n \times n}$ and $b \in \mathbb{R}^n$, we seek $x$ such that $Ax = b$. The Gaussian elimination process can be algorithmically implemented as follows:

def gaussian_elimination(A, b):
    n = len(A)
    # Forward elimination
    for k in range(n):
        for i in range(k+1, n):
            if A[k][k] == 0:
                raise ValueError("Zero pivot encountered.")
            alpha = A[i][k] / A[k][k]

            for j in range(k, n):
                A[i][j] -= alpha * A[k][j]
            b[i] -= alpha * b[k]

    # Backward substitution
    x = [0] * n
    for i in range(n-1, -1, -1):
        sum_ax = sum(A[i][j] * x[j] for j in range(i+1, n))
        x[i] = (b[i] - sum_ax) / A[i][i]

    return x

5.2.3 Example link

Consider the following system:

$$A = \begin{bmatrix} 1 & 4 & 7 \\ 2 & 5 & 8 \\ 3 & 6 & 11 \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}$$

We construct the elimination matrices $M_1$ and $M_2$ to systematically annihilate the sub-diagonal entries of the first and second columns, respectively:

$$M_1 = \begin{bmatrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -3 & 0 & 1 \end{bmatrix}, \quad M_2 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -2 & 1 \end{bmatrix}$$

Applying these transformations yields the upper triangular matrix $U$ and the updated right-hand side vector $\tilde{b}$:

$$U = M_2 M_1 A = \begin{bmatrix} 1 & 4 & 7 \\ 0 & -3 & -6 \\ 0 & 0 & 2 \end{bmatrix}, \quad \tilde{b} = M_2 M_1 b = \begin{bmatrix} 1 \\ -1 \\ 0 \end{bmatrix}$$

To solve the equivalent system $Ux = \tilde{b}$, we execute backward substitution:

$$x_3 = 0, \quad -3x_2 - 6(0) = -1 \implies x_2 = 1/3, \quad x_1 + 4(1/3) + 7(0) = 1 \implies x_1 = -1/3$$

Thus, the solution vector is exactly:

$$x = \begin{bmatrix} -1/3 \\ 1/3 \\ 0 \end{bmatrix}$$

5.2.4 LU Decomposition link

Assume that after $k-1$ steps, the intermediate matrix is partitioned as follows:

$$A^{(k-1)} = M_{k-1} \cdots M_1 A = \left[ \begin{array}{c|c} A_{11}^{(k-1)} & A_{12}^{(k-1)} \\ \hline 0 & A_{22}^{(k-1)} \end{array} \right] \begin{array}{l} k-1 \\ n-k+1 \end{array}$$

where $A_{11}^{(k-1)}$ is an upper triangular matrix, and $A_{22}^{(k-1)}$ is the remaining submatrix:

$$A_{22}^{(k-1)} = \begin{bmatrix} a_{kk}^{(k-1)} & \cdots & a_{kn}^{(k-1)} \\ \vdots & \ddots & \vdots \\ a_{nk}^{(k-1)} & \cdots & a_{nn}^{(k-1)} \end{bmatrix}$$

Assuming the pivot element $a_{kk}^{(k-1)} \neq 0$, we define the elimination matrix $M_k = I - \alpha^{(k)} e_k^T$, where the multiplier vector $\alpha^{(k)}$ is given by:

$$\alpha^{(k)} = \begin{bmatrix} 0 \\ \vdots \\ 0 \\ l_{k+1,k} \\ \vdots \\ l_{n,k} \end{bmatrix}, \quad \text{where } l_{ik} \triangleq \frac{a_{ik}^{(k-1)}}{a_{kk}^{(k-1)}}, \quad \text{for } k+1 \le i \le n$$

Applying $M_k$ advances the elimination process:

$$A^{(k)} = M_k A^{(k-1)} = \left[ \begin{array}{c|c} A_{11}^{(k)} & A_{12}^{(k)} \\ \hline 0 & A_{22}^{(k)} \end{array} \right] \begin{array}{l} k \\ n-k \end{array}$$

Finally, after $n-1$ steps, we obtain the upper triangular matrix $U$:

$$A^{(n-1)} = \left( \prod_{i=1}^{n-1} M_i \right) A = U$$

That is to say, the original matrix $A$ can be decomposed as:

$$A = \left( \prod_{i=1}^{n-1} M_i \right)^{-1} U = LU$$

Utilizing the property $(I - \alpha^{(k)} e_k^T)(I + \alpha^{(k)} e_k^T) = I$, we can explicitly construct $L$:

$$L = \prod_{i=1}^{n-1} (I + \alpha^{(i)} e_i^T) = I + \sum_{i=1}^{n-1} \alpha^{(i)} e_i^T$$

Here, $L$ is a unit lower triangular matrix (i.e., having $1$s on its diagonal part).

5.2.5 Algorithm link

import numpy as np

def lu_decomposition(A: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    # Performs LU decomposition of a square matrix A.
    # Returns L (unit lower triangular) and U (upper triangular).
    n = A.shape[0]
    L = np.eye(n) # Initialize L as identity matrix
    U = A.copy().astype(float) # Initialize U as a copy of A

    for k in range(n - 1):
        if U[k, k] == 0:
            raise ValueError(f"Zero pivot encountered at index {k}.")

        for i in range(k + 1, n):
            # Calculate the multiplier l_ik
            factor = U[i, k] / U[k, k]
            L[i, k] = factor

            # Apply the elementary row operation
            # U[i] = U[i] - factor * U[k]
            U[i, k:] -= factor * U[k, k:]

    return L, U

# --- Example Usage based on your notes ---
if __name__ == "__main__":
    A = np.array([[1, 4, 7],
                  [2, 5, 8],
                  [3, 6, 11]], dtype=float)
    L, U = lu_decomposition(A)

    print("Matrix L:\n", L)
    print("Matrix U:\n", U)
    # Verification: np.allclose(L @ U, A) will return True

5.3 Iterative Methods link

For solving a linear system $Ax = b$, we define the residual vector as $r \triangleq b - Ax$. There are two main types of iterative methods:

1. Stationary methods

2. Nonstationary methods

We can decompose the matrix $A$ as $A = D + L + U$, where:

• $D = \text{diag}(A)$

• $L = \text{strictly lower triangular part of } A$

• $U = \text{strictly upper triangular part of } A$

Thus, the system $Ax = b$ can be rewritten as $(D + L + U)x = b$.

5.3.1 Stopping Criteria for Iterative Methods link

Before introducing the specific methods, we define the standard stopping criteria for iterative algorithms. As shown in the general iteration loop, the process terminates when either the step size or the residual falls below predefined tolerances ($\varepsilon_1$ and $\varepsilon_2$):

$$\|x^{(k)} - x^{(k-1)}\| < \varepsilon_1 \quad \text{or} \quad \|b - Ax^{(k)}\| < \varepsilon_2$$

5.3.2 Jacobi Method link

Given an initial value $x^{(0)}$, the Jacobi method updates the solution iteratively:

$$D x^{(k)} = -(L + U)x^{(k-1)} + b, \quad k = 1, 2, \dots$$

Equivalently:

$$x^{(k)} = -D^{-1}(L + U)x^{(k-1)} + D^{-1}b$$

5.3.3 Algorithm link

import numpy as np

def jacobi_iteration(A: np.ndarray, b: np.ndarray, x0: np.ndarray,
                     tolerance: float = 1e-6, max_iter: int = 100) -> tuple[np.ndarray, int]:
    """
    Solves the linear system Ax = b using the Jacobi iterative method.
    """
    n = len(A)
    x = x0.copy().astype(float)
    x_new = np.zeros_like(x)

    for k in range(max_iter):
        for i in range(n):
            # Calculate the sum of a_ij * x_j for j != i
            s = sum(A[i, j] * x[j] for j in range(n) if j != i)
            # Update the i-th component of the new state
            x_new[i] = (b[i] - s) / A[i, i]

        # Calculate the residual norm to check for convergence
        if np.linalg.norm(x_new - x, ord=np.inf) < tolerance:
            return x_new, k + 1

        x = x_new.copy()

    print("Warning: Maximum iterations reached without convergence.")
    return x, max_iter

5.3.4 Gauss-Seidel Method link

Unlike the Jacobi method, which updates the solution vector only after computing all components, the Gauss-Seidel method utilizes the newly computed components $x_i^{(k)}$ as soon as they become available. Using the matrix decomposition $A = D + L + U$, the Gauss-Seidel iteration is formulated as:

$$(D + L)x^{(k)} = -Ux^{(k-1)} + b$$

By isolating $x^{(k)}$, we obtain the explicit iteration form:

$$x^{(k)} = \underbrace{-(D + L)^{-1}U}_{T} x^{(k-1)} + (D + L)^{-1}b$$

where $T$ serves as the iteration matrix.

5.3.5 Algorithm: Gauss-Seidel Iteration link

Input: Initial guess $x^{(0)}$, stopping tolerances $\varepsilon_1, \varepsilon_2$, maximum iterations $N$ For $k = 1 : N$ For $i = 1 : n$

$$x_i^{(k)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} \overbrace{x_j^{(k)}}^{\text{new}} - \sum_{j=i+1}^{n} a_{ij} \overbrace{x_j^{(k-1)}}^{\text{old}} \right)$$

End If $\|x^{(k)} - x^{(k-1)}\| < \varepsilon_1$ or $\|b - Ax^{(k)}\| < \varepsilon_2$ then Stop End End

5.3.6 Convergence of Fixed-Point Iterations link

Both the Jacobi and Gauss-Seidel methods are specific instances of fixed-point iterations, which can be generally expressed in the form:

$$x^{(k)} = Tx^{(k-1)} + c$$

To establish the conditions under which these sequences converge, we introduce the following foundational lemma and theorem regarding the iteration matrix $T$. Lemma. If the spectral radius $\rho(T)$ satisfies $\rho(T) < 1$, then the matrix $(I - T)$ is invertible (i.e., $(I - T)^{-1}$ exists), and its inverse can be expanded as a series:

$$(I - T)^{-1} = I + T + T^2 + \cdots = \sum_{i=0}^{\infty} T^i$$

Corollary. As a practical consequence, if $\|T\| < 1$ for any natural matrix norm, the iteration is mathematically guaranteed to converge.

5.3.7 Conjugate Gradient Method (CG Method) link

The Conjugate Gradient (CG) method is an advanced iterative algorithm designed to solve systems of linear equations $Ax = b$, where the matrix $A \in \mathbb{R}^{n \times n}$ is strictly symmetric and positive definite. A fundamental property of the CG method is that it constructs the solution space over a sequence of basis vectors $\{v^{(i)}\}_{i=1}^n$ which form an $A$-orthogonal (or conjugate) set. Mathematically, this conjugacy condition is expressed as:

$$\langle v^{(k)}, A v^{(i)} \rangle = 0, \quad \text{for } i = 1, 2, \dots, k-1$$

Algorithm: CG Method 0. Given an initial guess $x^{(0)}$ and a specified tolerance $\varepsilon$. 1. Initialize the residual and the initial search direction:

$$\begin{aligned} r^{(0)} &= b - A x^{(0)} \\ v^{(1)} &= r^{(0)} \end{aligned}$$

2. For $k = 1 : n$ Compute the step size $t_k$:

$$t_k = \frac{\langle r^{(k-1)}, r^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle}$$

Update the solution and the residual vector:

$$\begin{aligned} x^{(k)} &= x^{(k-1)} + t_k v^{(k)} \\ r^{(k)} &= r^{(k-1)} - t_k A v^{(k)} \end{aligned}$$

If $\|r^{(k)}\| < \varepsilon$, then stop. Compute the orthogonalization coefficient $s_k$:

$$s_k = \frac{\|r^{(k)}\|_2^2}{\|r^{(k-1)}\|_2^2}$$

Update the search direction for the next iteration:

$$v^{(k+1)} = r^{(k)} + s_k v^{(k)}$$

End

5.3.8 Projection Methods link

To generalize the principles underlying the CG method, we introduce projection methods. At each iteration $k$, we define two distinct subspaces: the search subspace $\mathcal{K}_k$ and the constraint subspace $\mathcal{L}_k$, with respective bases $V_k$ and $W_k$. We enforce the condition that the approximate solution resides in the affine space:

$$x^{(k)} \in x^{(0)} + \mathcal{K}_k \implies x^{(k)} = x^{(0)} + V_k y_k$$

A projection method where the constraint subspace is identical to the search subspace ($\mathcal{L}_k = \mathcal{K}_k$) is formally called the Galerkin method. For the CG method specifically, the search space is the Krylov subspace, defined as:

$$\mathcal{K}_k = \text{span}\{r^{(0)}, A r^{(0)}, A^2 r^{(0)}, \dots, A^{k-1} r^{(0)}\}$$

6 Approximation Theory link

While Krylov subspace methods are optimal for solving square deterministic systems, empirical observations often lead to overdetermined systems. Consider an experiment evaluating Hooke’s Law for a spring, modeled as:

$$e + kF = l$$

where $l$ is the measured length, $F$ is the applied force, and the constants $e$ (initial length) and $k$ (compliance) are unknown. Suppose an experiment yields the following dataset:

$F$	1	2	3	4	5
$l$	7.97	10.2	14.2	16.0	21.2

To determine the unknowns $e$ and $k$, we can formulate this as a linear system:

$$\begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \\ 1 & 5 \end{bmatrix} \begin{bmatrix} e \\ k \end{bmatrix} = \begin{bmatrix} 7.97 \\ 10.2 \\ 14.2 \\ 16.0 \\ 21.2 \end{bmatrix}$$

Because this system has more equations than unknowns, an exact solution rarely exists. This naturally leads to the necessity of Discrete Least Squares Approximation. Assume we have a general group of data points $\{(x_i, y_i)\}_{i=1}^n$. We wish to represent these data using a predefined functional form. If we assume a linear function $y = ax + b$, the objective shifts to answering a fundamental question: how do we determine the optimal parameters $a$ and $b$ to minimize the discrepancy between the model and the data?

6.1 Discrete Least Squares Approximation link

6.1.1 Error Definitions in Approximation link

To quantify the discrepancy between the measured data points $\{(x_i, y_i)\}_{i=1}^n$ and our proposed linear model $y = ax + b$, we define three fundamental types of measurement errors:

1. Maximum Error (Chebyshev norm): $E_\infty(a,b) \triangleq \max_{1 \le i \le n} |y_i - (ax_i + b)|$

2. Absolute Error (L1 norm): $E_1(a,b) \triangleq \sum_{i=1}^n |y_i - (ax_i + b)|$

3. Least Squares Error (L2 norm): $E_2(a,b) \triangleq \left( \sum_{i=1}^n [y_i - (ax_i + b)]^2 \right)^{1/2}$

In general, we determine the optimal parameters $a$ and $b$ based on minimizing the smallest error. The most analytically tractable choice is the Best Least Squares approximation. To find the minimum, we set the partial derivatives of the squared error to zero:

$$\frac{\partial (E_2)^2}{\partial a} = 0, \quad \frac{\partial (E_2)^2}{\partial b} = 0$$

This differentiation yields a system of equations:

$$\begin{aligned} -2 \sum_{i=1}^n \big(y_i - (ax_i + b)\big)x_i = 0 &\implies a\sum_{i=1}^n x_i^2 + b\sum_{i=1}^n x_i = \sum_{i=1}^n x_i y_i \\ -2 \sum_{i=1}^n \big(y_i - (ax_i + b)\big) = 0 &\implies a\sum_{i=1}^n x_i + bn = \sum_{i=1}^n y_i \end{aligned}$$

Solving this $2 \times 2$ linear system provides the exact analytical expressions for $a$ and $b$.

6.1.2 The Normal Equations link

To answer the question of determining the optimal parameters $x = \begin{bmatrix} b \\ a \end{bmatrix}$ (or $x = \begin{bmatrix} e \\ k \end{bmatrix}$ in our Hooke’s Law example) for an overdetermined system $Ax = y$, we must acknowledge that an exact solution rarely exists. Instead, we seek to minimize the $L_2$ norm of the residual vector, $\|Ax - y\|_2^2$. Geometrically, the optimal approximation occurs when the residual vector $(y - Ax)$ is orthogonal to the column space of $A$. This orthogonality principle is algebraically enforced by multiplying both sides of our system by the transpose of the matrix, $A^T$:

$$A^T A x = A^T y$$

These are known as the Normal Equations. Provided that the columns of $A$ are linearly independent (which is typical for distinct data points in experimental models), the square matrix $A^T A$ is symmetric and positive definite. This guarantees that its inverse exists, leading to the unique least squares solution:

$$x = (A^T A)^{-1} A^T y$$

This derivation conceptually unites our previous discussions on projection methods, orthogonalization, and solving systems involving symmetric positive definite matrices.

6.1.3 Polynomial Least Squares Approximation link

If we consider approximating the function using a polynomial of degree $m$ for the data $\{(x_i, y_i)\}_{i=1}^N$ (where $0 < m < N-1$):

$$P_m(x) = a_m x^m + a_{m-1} x^{m-1} + \cdots + a_1 x + a_0$$

The total least squares error becomes:

$$E_2(a_0, \dots, a_m) \triangleq \sum_{i=1}^N \big(y_i - P_m(x_i)\big)^2$$

Similarly, we set $\frac{\partial E_2}{\partial a_k} = 0$ for $k = 0, 1, \dots, m$ to find the minimum value. This leads to a linear system of equations $A u = b$:

$$\begin{bmatrix} N & \sum x_i & \cdots & \sum x_i^m \\ \sum x_i & \sum x_i^2 & \cdots & \sum x_i^{m+1} \\ \vdots & \vdots & \ddots & \vdots \\ \sum x_i^m & \sum x_i^{m+1} & \cdots & \sum x_i^{2m} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_m \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum y_i x_i \\ \vdots \\ \sum y_i x_i^m \end{bmatrix}$$

Note: You can use iterative methods like the CG method to solve this linear system. In software environments like MATLAB, built-in functions such as polyfit or lsqlin can be used to efficiently approximate these linear functions.

6.2 Continuous Least Squares Approximation link

Suppose $f \in C[a,b]$. We want to approximate it through a polynomial of degree at most $n$, $P_n(x) = \sum_{k=0}^n a_k x^k$. The error function is defined continuously as:

$$E_2(a_0, a_1, \dots, a_n) \triangleq \int_a^b [f(x) - P_n(x)]^2 dx = \|f - P_n\|_2^2$$

To minimize this error, for each index $j = 0, 1, \dots, n$, we set $\frac{\partial E_2}{\partial a_j} = 0$:

$$\begin{aligned} \frac{\partial}{\partial a_j} \int_a^b \big(f^2 - 2f P_n + P_n^2\big) dx &= -2 \int_a^b x^j f(x) dx + 2 \sum_{k=0}^n a_k \int_a^b x^{j+k} dx = 0 \\ \implies \sum_{k=0}^n a_k \int_a^b x^{j+k} dx &= \int_a^b x^j f(x) dx \end{aligned}$$

This leads to the continuous counterpart of our matrix system $A u = b$:

$$\begin{bmatrix} \int_a^b 1 \,dx & \int_a^b x \,dx & \cdots & \int_a^b x^n \,dx \\ \int_a^b x \,dx & \int_a^b x^2 \,dx & \cdots & \int_a^b x^{n+1} \,dx \\ \vdots & \vdots & \ddots & \vdots \\ \int_a^b x^n \,dx & \int_a^b x^{n+1} \,dx & \cdots & \int_a^b x^{2n} \,dx \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} \int_a^b f(x) \,dx \\ \int_a^b x f(x) \,dx \\ \vdots \\ \int_a^b x^n f(x) \,dx \end{bmatrix}$$

Remark: The monomials $\{1, x, x^2, \dots, x^n\}$ act as the basis for the subspace of $P_n(x)$. However, this basis can lead to ill-conditioned matrices (the Hilbert matrix). We can consider another subspace with orthogonal bases $\{\phi_i(x)\}_{i=1}^n$ for better numerical stability. Let $\Pi_n$ be the set of all polynomials of degree at most $n$. Any polynomial $P \in \Pi_n$ can be represented as a linear combination of these basis functions: $P = \sum_{i=1}^n c_i \phi_i(x)$.

6.3 Orthogonal Polynomials link

6.3.1 Legendre Polynomials link

A simple basis function set for approximation on the interval $[a,b] = [-1, 1]$ is the monic Legendre polynomials:

$$\begin{aligned} P_0(x) &= 1 \\ P_1(x) &= x \\ P_2(x) &= x^2 - \frac{1}{3} \\ P_3(x) &= x^3 - \frac{3}{5}x \\ P_4(x) &= x^4 - \frac{6}{7}x^2 + \frac{3}{35} \end{aligned}$$

These satisfy the orthogonality condition: $\int_{-1}^1 P_k(x) P_j(x) dx = 0$ if $k \neq j$.

6.3.2 Chebyshev Polynomials link

Defined as $T_n(x) \triangleq \cos(n \cos^{-1} x)$ for $n \ge 0$ and $x \in [-1, 1]$. The sequence $\{T_n(x)\}_{n=1}^N$ is orthogonal with respect to the weighting function $w(x) = \frac{1}{\sqrt{1-x^2}}$:

$$\int_{-1}^1 T_m(x) T_n(x) w(x) dx = 0 \quad \text{if } m \neq n$$

Let $\theta = \cos^{-1} x \implies x = \cos \theta$. Then $T_n(\theta) = \cos(n\theta)$. Utilizing trigonometric identities:

$$\begin{aligned} T_{n+1}(\theta) &= \cos((n+1)\theta) = \cos(n\theta)\cos\theta - \sin(n\theta)\sin\theta \\ T_{n-1}(\theta) &= \cos((n-1)\theta) = \cos(n\theta)\cos\theta + \sin(n\theta)\sin\theta \end{aligned}$$

Adding these yields the recurrence relation: $T_{n+1}(\theta) + T_{n-1}(\theta) = 2\cos(n\theta)\cos\theta$. Converting back to $x$:

$$T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x) \quad \text{for } n \ge 1$$

Generating the first few polynomials:

$$\begin{aligned} T_0(x) &= 1 \\ T_1(x) &= x \\ T_2(x) &= 2x T_1(x) - T_0(x) = 2x^2 - 1 \\ T_3(x) &= 2x T_2(x) - T_1(x) = 4x^3 - 3x \end{aligned}$$

7 Approximating Eigenvalues link

7.1 Singular Value Decomposition (SVD) link

First, we verify that $A_k$ achieves this error bound:

$$\|A - A_k\|_2 = \left\| U(A - A_k)V^T \right\|_2 = \|\text{diag}(0, \dots, 0, \sigma_{k+1}, \sigma_{k+2}, \dots, \sigma_p)\|_2 = \sigma_{k+1}$$

Since $\text{rank}(\hat{A}) = k$, for all $\hat{A} \in \mathbb{R}^{m \times n}$, there exist orthonormal vectors $\{x_1, \dots, x_{n-k}\}$ such that $\text{Null}(\hat{A}) = \text{span}\{x_1, \dots, x_{n-k}\}$. (If $y \in \text{Null}(\hat{A})$, $y = \sum_{i=1}^{n-k} c_i x_i$, then $\hat{A}y = 0$). Assume we construct the intersection space $Z = \text{span}\{x_1, x_2, \dots, x_{n-k}\} \cap \text{span}\{v_1, v_2, \dots, v_{k+1}\}$. Because the sum of the dimensions is $(n-k) + (k+1) = n+1 > n$, this intersection space is non-empty ($Z \neq \{0\}$). Then $\exists z \in Z$ with $\|z\|_2 = 1$ such that $\hat{A}z = 0$. We evaluate $Az$:

$$Az = \left( \sum_{i=1}^{k+1} \sigma_i u_i v_i^T \right) z = \sum_{i=1}^{k+1} \sigma_i (v_i^T z) u_i$$

$$\implies \|Az\|_2^2 = \sum_{i=1}^{k+1} \sigma_i^2 (v_i^T z)^2$$

Therefore, the error is bounded from below:

$$\|A - \hat{A}\|_2^2 \ge \|(A - \hat{A})z\|_2^2 = \|Az\|_2^2 = \sum_{i=1}^{k+1} \sigma_i^2 (v_i^T z)^2 \ge \sigma_{k+1}^2 \sum_{i=1}^{k+1} (v_i^T z)^2 = \sigma_{k+1}^2 \|z\|_2^2 = \sigma_{k+1}^2$$

Hence, $\|A - \hat{A}\|_2 \ge \sigma_{k+1}$.

Remark: If $\text{rank}(\hat{A}) = r_\varepsilon$ with a threshold $\varepsilon = M \|A\|_2$, then $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_{r_\varepsilon} > \varepsilon > \sigma_{r_\varepsilon+1} \ge \dots \ge \sigma_p \ge 0$. We say $\hat{A} \sim A$ within the significance of error $\varepsilon$. Geometric Interpretation of Singular Values (Hyper-space): Consider $A = \begin{bmatrix} 0.96 & 1.72 \\ 1.72 & 0.96 \end{bmatrix}$. The SVD is $A = U \Sigma V^T = \begin{bmatrix} \frac{3}{5} & -\frac{4}{5} \\ \frac{4}{5} & \frac{3}{5} \end{bmatrix} \begin{bmatrix} 3 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \frac{4}{5} & \frac{3}{5} \\ \frac{3}{5} & -\frac{4}{5} \end{bmatrix}^T$. If we define a set $E = \{y \mid y = Ax, \forall x \in \mathbb{R}^2, \|x\|_2 = 1\}$, the singular values represent the lengths of the semi-axes of the resulting hyper-ellipse.

7.2 Orthogonal Projection link

Let $S = \text{span}\{v\}$ for a vector $v \in \mathbb{R}^n$. We seek a matrix $P$ such that for any $x \in \mathbb{R}^n$, $Px \in S$ and the residual $(I - P)x \in S^\perp$. Let $P \triangleq \frac{v v^T}{\|v\|_2^2} = \left(\frac{v}{\|v\|_2}\right) \left(\frac{v}{\|v\|_2}\right)^T$. Therefore, for any $x \in \mathbb{R}^n$:

$$Px = v \left( \frac{v^T x}{\|v\|_2^2} \right)$$

This matrix $P$ satisfies:

$$\begin{cases} P^2 = P \cdot P = P \\ P^T = P \end{cases}$$

so we call $P$ an orthogonal projection. If $V = [v_1, v_2, \dots, v_k]$ with orthonormal basis $v_i$ for $S$, then $P = V V^T$ is the unique orthogonal projection onto $S$. Projections from SVD: Let $A = [U_r, \bar{U}_r] \Sigma [V_r, \bar{V}_r]^T$ be the SVD of $A$, where $\text{rank}(A) = r$. $U_r = [u_1, \dots, u_r]$, $\bar{U}_r = [u_{r+1}, \dots, u_m]$, $\Sigma = \text{diag}(\sigma_1, \dots, \sigma_p)$, $p = \min(m,n)$. $V_r = [v_1, \dots, v_r]$. Then we have the following results:

1. $P = V_r V_r^T$ is the orthogonal projection onto $N(A)^\perp = R(A^T)$.

2. $P = \bar{V}_r \bar{V}_r^T$ is the orthogonal projection onto $N(A)$.

3. $P = U_r U_r^T$ is the orthogonal projection onto $R(A)$.

4. $P = \bar{U}_r \bar{U}_r^T$ is the orthogonal projection onto $R(A)^\perp = N(A^T)$.

7.2.1 Application: Distance Between Subspaces link

To measure the distance between two subspaces $S_1, S_2 \in \mathbb{R}^n$ with $\dim(S_1) = \dim(S_2)$. Assume $P_1, P_2$ are two orthogonal projections for $S_1, S_2$ respectively. Then:

$$\text{dist}(S_1, S_2) = \|P_1 - P_2\|_2$$

Example: Let $S_1 = \text{span}\{x\} \in \mathbb{R}^2$ where $x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$, and $S_2 = \text{span}\{y\} \in \mathbb{R}^2$ where $y = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}$. Assume $\|x\|_2 = \|y\|_2 = 1$. The angle between them is defined by $x^T y = \cos\theta$, for $0 \le \theta \le \frac{\pi}{2}$.

$$P_1 = \frac{x x^T}{x^T x} = x x^T, \quad P_2 = \frac{y y^T}{y^T y} = y y^T$$

The difference is $P_1 - P_2 = x x^T - y y^T$. Through SVD analysis, this difference matrix has non-zero singular values equal to $\sin\theta$.

$$\therefore \text{dist}(S_1, S_2) = \|P_1 - P_2\|_2 = \sin\theta$$

Solving Non-linear Vector Functions Consider solving a non-linear vector function $F(x) = 0$, where $F: \mathbb{R}^n \to \mathbb{R}^m$. For example, let $F: \mathbb{R}^n \to \mathbb{R}^1$ be a quadratic function:

$$F(x) = \frac{1}{2} x^T A x - x^T b, \quad A \in \mathbb{R}^{n \times n}, \ x, b \in \mathbb{R}^n$$

Using the gradient method (such as CG) for solving $Ax = b$, we essentially find the roots of the derivative:

$$\frac{\partial F}{\partial x} = Ax - b = 0$$

to ensure the minimum of $F(x)$. To generalize this to broader non-linear systems, we must establish rigorous definitions for multivariable derivatives.

8 Multivariable Calculus link

8.1 F-Differentiable (Fréchet Derivative) link

Recall Newton’s method in 1D: for $f(x) = 0$, the update rule is $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$. To extend this, let $D \subset \mathbb{R}^n$ be a given domain and let $F: D \subset \mathbb{R}^n \to \mathbb{R}^m$ be a given non-linear function.

Remark:

1. If $A$ exists, $A$ is the Jacobian matrix. $A = (a_{ij})$ for $1 \le i \le m$, $1 \le j \le n$, where $a_{ij} = \frac{\partial F_i}{\partial x_j}$ and $F = [F_1, F_2, \dots, F_m]^T$.

2. We say $F$ is continuously differentiable, denoted by $F \in C^1(D)$, if $F'(x)$ is continuous.

Example 1:

$$f(x_1, x_2) = \begin{cases} x_1 & \text{if } x_2 = 0 \\ x_2 & \text{if } x_1 = 0 \\ 1 & \text{otherwise} \end{cases}$$

where $f: \mathbb{R}^2 \to \mathbb{R}^1$. The partial derivatives $\frac{\partial f}{\partial x_1}(0)$ and $\frac{\partial f}{\partial x_2}(0)$ exist. However, $f$ is not F-differentiable at $0$ because it is not continuous there.

8.2 G-Differentiable (Gâteaux Derivative) link

Example 2:

$$f(x_1, x_2) = \begin{cases} 0 & \text{if } x_1 = 0 \\ \frac{2x_2 e^{-x_1^2}}{x_2^2 + 2^{-2x_1^2}} & \text{if } x_1 \neq 0 \end{cases}$$

$f$ has a G-derivative at $0$, but $f$ is discontinuous at $x = 0$. Example 3:

$$f(x_1, x_2) = \begin{cases} 0 & \text{if } (x_1, x_2) = (0,0) \\ \frac{x_2 (x_1^2 + x_2^2)^{3/2}}{x_1^4 + x_2^2} & \text{if } (x_1, x_2) \neq (0,0) \end{cases}$$

This is a classic function used to show that $f$ has a G-derivative at $x=0$, but not an F-derivative at $x=0$.

8.2.1 Exercise link

• Lemma 1: If $F: D \subset \mathbb{R}^n \to \mathbb{R}^m$ is F-differentiable, then $F$ is continuous.

• Lemma 2: If $F$ is F-differentiable at $x \in D$, then $\frac{\partial F_i}{\partial x_j}$ exists at $x \in D$ and forms the Jacobian matrix.

9 Systems of Nonlinear Equations link

Solving a nonlinear vector function $F(x) = 0$, where $F: \mathbb{R}^n \to \mathbb{R}^n$, requires extrapolating iterative scalar methods into multi-dimensional spaces.

9.1 Fixed-Point Iteration for Systems link

Let $G: D \subset \mathbb{R}^n \to \mathbb{R}^n$ where $D$ is an open set. The iterative process identifies the fixed point $x^* \in D$ such that $x^* = G(x^*)$ by continuously applying the recurrence relation:

$$x_{k+1} = G(x_k)$$

A fixed point $x^*$ is defined theoretically as a “point of attraction” if there exists a bounded continuous region $B \subset D$ encompassing $x^*$ such that $G(B) \subset B$ and any initialization $x_0 \in B$ results in $x_k \to x^*$. If $G \in C^1(D)$ and the spectral radius of $G'(x^*)$ is strictly less than 1, the iteration converges reliably. This is intrinsically linked to the Contraction Mapping Theorem, stating that if $\|G(x) - G(y)\| \le \alpha \|x - y\|$ with $0 < \alpha < 1$, then the sequence will be a Cauchy sequence and inevitably converge linearly.

9.2 Newton’s Method for Systems link

Consider $F: D \subset \mathbb{R}^n \to \mathbb{R}^n$ where $F(x^*) = 0$. Using the multidimensional Taylor expansion near $x_k$:

$$F(x) = F(x_k) + F'(x_k)(x - x_k) + R(x, x_k)$$

Setting the linear tangent approximation $T_k(x_{k+1}) = 0$ provides the iterative update rule:

$$0 = F(x_k) + F'(x_k)(x_{k+1} - x_k) \implies x_{k+1} = x_k - [F'(x_k)]^{-1} F(x_k)$$

In practical computational implementations, explicitly calculating the matrix inverse is severely inefficient. Instead, we structurally define a step vector $s_k$ and solve the equivalent linear system:

$$F'(x_k) s_k = -F(x_k)$$

$$x_{k+1} = x_k + s_k$$

This method achieves local quadratic convergence provided the Jacobian $F'(x^*)$ is definitively nonsingular. To eliminate the significant overhead of computing analytic derivatives, the Jacobian matrix can be programmatically approximated utilizing forward finite differences:

$$A_{ij}(x) \triangleq \frac{f_i(x + h e_j) - f_i(x)}{h}$$

where $e_j$ is the canonical standard basis vector and the perturbation factor $h$ is generally scaled to $h \approx \sqrt{\epsilon_{machine}}$ to perfectly balance theoretical truncation error against numerical round-off error.

9.3 Inexact Newton’s Method and GMRES link

Solving the exact linear system $F'(x_k) s_k = -F(x_k)$ incurs an extreme $\mathcal{O}(n^3)$ cost per iteration. The Inexact Newton’s Method deliberately relaxes this strict condition, solving the system approximately such that the residual decreases up to a designated tolerance limit:

$$\|F'(x_k) s_k + F(x_k)\| \le \eta_k \|F(x_k)\|$$

For highly asymmetric and sparse Jacobian configurations, the Generalized Minimum Residual (GMRES) method becomes strictly necessary. GMRES approximates the iterative solution by effectively minimizing the residual norm over a continuously expanding Krylov subspace $\mathcal{K}_k = \text{span}\{r_0, A r_0, \dots, A^{k-1} r_0\}$. The Finite Difference GMRES (FDGMRES) variant operates entirely matrix-free by directly estimating the Jacobian-vector products:

$$F'(x)v \approx \frac{F(x + hv) - F(x)}{h}$$

FDGMRES methodically constructs an orthonormal basis $V_k$ using the Arnoldi process, converting the high-dimensional problem into minimizing $\|\beta e_1 - H_k y_k\|$ across a much smaller Upper-Hessenberg matrix $H_k$.

10 Numerical Solutions for Boundary Value Problems link

To accurately model physical displacement distributions, we utilize specialized discretization frameworks.

10.1 Finite Difference Method (FDM) link

To numerically solve a boundary value problem, the Finite Difference Method discretizes the domain into a uniform grid. Consider a theoretical physical string displacement problem constrained by:

$$-u''(x) = f(x), \quad 0 < x < 1$$

$$u(0) = u(1) = 0$$

We partition the interval uniformly, defining $x_i = i h$ for $i = 0, 1, \dots, n$ where $h = \frac{1}{n}$. The continuous second derivative is analytically approximated employing the central difference scheme:

$$u''(x_i) = \frac{u(x_{i-1}) - 2u(x_i) + u(x_{i+1})}{h^2} + \mathcal{O}(h^2)$$

Substituting this localized approximation directly into the governing differential equation yields:

$$-\left( \frac{u_{i-1} - 2u_i + u_{i+1}}{h^2} \right) = f(x_i)$$

This explicitly translates the continuum calculus problem into an algebraically structured linear system $Au = b$. Because each node connects exclusively to its immediate adjacent neighbors, the resulting interaction matrix $A$ becomes a definitively tridiagonal matrix scaled by $\frac{1}{h^2}$, with diagonal elements measuring $2$ and sub-diagonal elements evaluating to $-1$.

10.2 Finite Element Method (FEM) link

The Finite Element Method approaches the differential equation problem by reformulating it integrally into its “weak form”. Multiplying the initial equation by a generic test function $v(x)$ and integrating continuously over the entire domain generates:

$$-\int_0^1 u'' v \, dx = \int_0^1 f v \, dx$$

Performing integration by parts directly on the left side transfers one structural derivative onto the specific test function. Given the inherent boundaries $v(0)=v(1)=0$, the boundary integral terms identically evaluate to zero, concluding with:

$$\int_0^1 u' v' \, dx = \int_0^1 f v \, dx$$

Instead of mathematically solving for isolated discrete points, the global solution space is systematically projected onto a highly localized finite-dimensional subspace formed by piecewise linear basis “hat” functions $\phi_j(x)$. Letting $u(x) \approx \sum_{j=1}^{n-1} c_j \phi_j(x)$ and choosing $v(x) = \phi_i(x)$, we assemble a discrete linear system $Ac = b$:

$$A_{ij} = \int_0^1 \phi_i'(x) \phi_j'(x) \, dx, \quad b_i = \int_0^1 f(x) \phi_i(x) \, dx$$

Because $A$ operates as structurally symmetric and definitively positive definite, it perfectly pairs with Conjugate Gradient solvers.