【笔记】Numerical Computation

小毛龙

老师写的ppt太狗屎了
我自己整理一下笔记好了

小毛龙

Lecture 01. Introduction

The Classification of Errors

Modeling Error
Measurement Error
Truncation Error
Round-off Error

Criteria of Evaluating Computation Algorithms

Complexity (time, space)
Accuracy
Stability

小毛龙

Lecture 02. Linear Systems

Solution of Linear Systems

A linear system can have:
- no solution (inconsistent)
- exact one solution (consistent)
- infinite solutions (consistent)
The basic method for solving a system of linear equations is to replace the given system by a new system that has the same solution set but which is easier to solve.
Since the rows of an augmented matrix correspond to the equations in the associated system. new systems is generally obtained in a series of steps by applying the following three types of operations to eliminate unknowns systematically. These are called elementary row operations:
1. Multiply an equation through by an nonzero constant
2. Interchange two equations
3. Add a constant multiple of one equation to another

Algorithms to Solve Linear Systems (Direct Methods)

Gauss-Jordan Elimination

I don't see the necessity to write the detailed procedure here.

Cramer's Rule

Find Inverse (A^{-1})

It is easy to consider that the solution to Ax = b can be written as x = A^{-1}b.
The problem is, how do we obtain the inverse A^{-1} ?

Simple Way

\begin{aligned} A &= \begin{bmatrix} a_{11} &\cdots &a_{1n}\\ \vdots &\vdots &\vdots\\ a_{m1} &\cdots &a_{mn}\\ \end{bmatrix}\\ &\longrightarrow \begin{bmatrix} a_{11} &\cdots &a_{1n} &| &1 &\cdots &1\\ \vdots &\vdots &\vdots &| &\vdots &\vdots &\vdots\\ a_{m1} &\cdots &a_{mn} &| &1 &\cdots &1\\ \end{bmatrix} \longrightarrow \begin{bmatrix} & &| & & &\\ &A &| & &I &\\ & &| & & &\\ \end{bmatrix}\\ &\longrightarrow \begin{bmatrix} & &| & & &\\ &I &| & &A^{-1} &\\ & &| & & &\\ \end{bmatrix} \end{aligned}

LU Decomposition

If all n leading principal minors of the n\times n matrix A are non-singular, then A has exact one LU-decomposition:

A = LU

Where L is a lower triangular matrix and U is an upper triangular matrix.

Algorithms to Solve Linear Systems (Iterative Methods)

\begin{aligned} A &= D - L -U, \begin{cases} D &= \{a_{ii} | a_{ii} \in A, i \in [1, n]\} \quad \text{(the diagonal line)}\\ L &= \{a_{ij} | a_{ij} \in A, i > j,i \in [1, n]\} \quad \text{(the neg. of lower part)}\\ U &= \{a_{ij} | a_{ij} \in A, i < j,i \in [1, n]\} \quad \text{(the neg. of upper part)} \end{cases}\\ \end{aligned}

Jacobi Method

T = D^{-1}(L+U), \quad c= (D-L)^{-1}b

Gauss-Seidel Method

T = (D-L)^{-1}U, \quad c= (D-L)^{-1}b

x^{(k)} = Tx^{(k-1)} + c, \text{where } k \in \N^+

Vector and Matrix Norms

Vector Norms

On a vector space V, a norm is a function ||v|| from the space V to \R^{+} that obeys these three postulates:

v \ne 0 \rightarrow ||v|| > 0
||\lambda v|| = |\lambda|\times ||v||, \lambda \in \R
||x + y|| \le ||x|| + ||y||

Then,
||x||_{l_p} = \bigg(\sum_{i=1}^{n} |x_i|^p\bigg)^{1/p}

Thus,
\begin{aligned} ||x||_1 &= \sum_{i = 1}^{n} |x_i|\\ ||x||_2 &= \sqrt{\sum_{i = 1}^{n} |x_i|^2}\\ ||x||_\infty &= \max_{i \in [1, n]} |x_i|\\ \end{aligned}

Matrix Norms

\begin{aligned} ||A||_1 &= \max_{j \in [1, m]}\sum_{i = 1}^{n} |a_{ij}| \quad (\text{largest column abs. sum})\\ ||x||_2 &= \sqrt{\ \lambda_{max}(A^{T}A)}\\ ||x||_\infty &= \max_{i \in [1, n]}\sum_{j = 1}^{m} |a_{ij}| \quad (\text{largest row abs. sum}) \end{aligned}

小毛龙

Lecture 03. Nonlinear Systems

Algorithms to Solve Nonlinear Equation

Bisection Method

Base on the fact if f(a)\times f(b) \le 0, and the function is continuous, then there must be a c that f(c) = 0, c \in [a, b].

We just repeatedly test if c lies in [a, \frac{a+b}2] or [\frac{a+b}2, b]. This algorithm is guaranteed to converge to the root.

However, this method is limited that it cannot solve the case where a function has a tangent of x = 0 at somewhere.

Newton's Method

x_{k+1} = x_{k}-\frac{f(x_k)}{f'(x_k)}, \ k\in \N^+

where x_0 is your choice (somewhere reasonably close to the root).

Conditions for Netown's Method to converge:
f \in C^2[a, b]

f(a)f(b) < 0

\forall x_i, x_j \in [a, b], f'(x_i) \ne 0, \text{sign}(f''(x_i)) = \text{sign}(f''(x_j))

\exist x_k \in [a, b], f(x_k)f''(x_k) > 0

Secant Method

This method is usually used when f'(x) is hard to calculate.

We approximate f'(x) by

f'(x_n) \approx \frac{f(x_n)-f(x_{n-1})}{x_n - x_{n-1}}

Then by Newton's Method:

x_{k+1} = x_{k}-\frac{f(x_k)}{f'(x_k)} = x_k -f(x_k)\times \bigg(\frac{f(x_k)-f(x_{k-1})}{x_k - x_{k-1}}\bigg)^{-1}

Descent Method

Descent Method is a slight variation of Newton's Method.
The basic principle is that if x_{k+1} calculated from x_k that |f(x_{k+1})| \ge |f(x_k)|, we try to find a x_j \in (x_k, x_{k+1}) that |f(x_j)| < |f(x_k)|.

x_p = \lambda(x_k - \frac{f(x_k)}{f'(x_k)}) + (1-\lambda) x_k = x_k - \lambda\frac{f(x_k)}{f'(x_k)}

\lambda \in (0, 1)

If \lambda = 1, then it is the original Newton's Method.

Steffensen's Method

A slight variation of Secant Method.
This method approximate f'(x) by
f'(x_n) \approx \frac{f(x_n + f(x_n ))-f(x_{n-1})}{f(x_n)}
Thus
x_{k+1} = x_{k}-\frac{f(x_k)}{f'(x_k)} = x_k -f(x_k)\times \bigg(\frac{f(x_k + f(x_k))-f(x_{k-1})}{f(x_k)}\bigg)^{-1}

Algorithms to Solve Nonlinear Systems

Consider a system of nonlinear equations:
\begin{aligned} &f_1(x_1, x_2,\dots, x_n) &= 0\\ &f_n(x_1, x_2,\dots, x_n) &= 0\\ &\vdots \\ &f_n(x_1, x_2,\dots, x_n) &= 0\\ \end{aligned}

Fixed Points of Nonlinear Equation

小毛龙

Lecture 04. Approximating Functions

Weierstrass Theorem

If f(x) is a continuous function on [a, b] and if \varepsilon > 0, then there is a polynomial p(x) satisfying
|f(x)-p(x)| \le \varepsilon
on the interval [a, b].

Lagrange Interpolation Polynomial

If a = x_0 < x_1 < \cdots < x_n = b are distinct real numbers on [a, b], then for arbitrary values y_0, y_1, \dots, y_n, there is a unique polynomial L_n(x) of degree n such that
L_n(x_i) = y_i, \quad i \in \N
We express L_n(x) as
L_n(x) = \sum_{k=0}^n y_kl_k(x)
where l_k(x) of degree n is
l_k(x) = \prod_{i = 0}^n \frac{x-x_i}{x_k-x_i} = \frac{(x-x_0)(x-x_1)\cdots(x-x_n)}{(x_k-x_0)(x_k-x_1)\cdots(x_k-x_n)}, l_i(x_j) = \begin {cases} 1, \quad i = j\\ \\ 0, \quad i \ne j \end {cases}

Example
x_0 = 1, x_1 = 2, x_2 = 4
y_0 = 2, y_1 = 1, y_2 = 3
\begin{aligned} L_2(x) &= y_0l_0(x) + y_1l_1(x) + y_2l_2(x)\\ &= 2\times \frac{(x-x_1)(x-x_2)}{(x_0 - x_1)(x_0-x_2)} + 1\times \frac{(x-x_0)(x-x_2)}{(x_1 - x_0)(x_1-x_2)} + 3\times \frac{(x-x_0)(x-x_1)}{(x_2 - x_0)(x_2-x_1)}\\ &= \frac 2 3 x^2-3x+\frac 1 3 \end{aligned}

Bernstein Polynomial

Spline Methods of Interpolation

Polynomials are the most common choice of interpolants because they are easy to:

Evaluate
Differentiate
Integrate

Linear Interpolation

Quadratic Interpolation

小毛龙

Lecture 05. Numerical Differentiation and Integration

Numerical Differentiation

Forward Difference

\begin{aligned} f'(x) = \lim_{h \to 0} \frac{f(x+h)-f(x)}{h} \end{aligned}

Backward Difference

\begin{aligned} f'(x) = \lim_{h\to 0} \frac{f(x) - f(x-h)}{h} \end{aligned}

Central Difference

\begin{aligned} f'(x) &= \lim_{h\to 0} \frac{1}{2} \bigg(\frac{f(x+h)-f(x)}{h} + \frac{f(x) - f(x-h)}{h}\bigg)\\ &= \lim_{h\to 0} \frac{f(x+h) - f(x-h)}{2h} \end{aligned}

First Three-Point Formula

\begin{aligned} f'(x) = \lim_{h\to 0} \frac{-3f(x) - 4f(x+h)- f(x+2h)}{2h} \end{aligned}

Second Three-Point Formula

\begin{aligned} f'(x) = \lim_{h\to 0} \frac{f(x+h)-f(x-h)}{2h} \end{aligned}

Numerical Integration

Rectangular Rule

\begin{aligned} \int_{x_0}^{x_0+h} f(x)\text d x\approx h\times f\bigg(\frac{x_0 + (x_0+h)}{2}\bigg) \end{aligned}

Trapzoid Rule

This method treats each subinterval as a trapzoid as bottom and top be f(x_0) and f(x_0+h) with a height of h.

\begin{aligned} \int_{x_0}^{x_0+h} f(x)\text d x\approx \frac{h}{2}(f(x_0)+f(x_0+h) \end{aligned}

Simpson's Rule

\begin{aligned} \int_{x_0}^{x_0+2h} f(x)\text d x\approx \frac{h}{3}\bigg[f(x_0) + 4f(x_0+h)+f(x_0+2h)\bigg] \end{aligned}

Composite Rectangular Rule

\begin{aligned} \int_{x_0}^{x_1} f(x)\text d x &\approx \sum_{k = 0}^{n-1} h\times f\bigg(\frac{(x_k) + (x_{k+1})}{2}\bigg)\\ h &= \frac{x_1 - x_0}{n}, x_k = x_0 + kh \end{aligned}

Composite Trapzoid Rule

\begin{aligned} \int_{x_0}^{x_1} f(x)\text d x &\approx \frac{h}{2}\times \bigg[f(x_0) + 2\sum_{k = 1}^{n-1} f(x_k) + f(x_1)\bigg]\\ h &= \frac{x_1 - x_0}{n}, x_k = x_0 + kh \end{aligned}

Composite Simpson's Rule

This method requires n must be even.
\begin{aligned} \int_{x_0}^{x_0+2h} f(x)\text d x &\approx \frac{h}{3}\bigg[f(x_0) + 4\sum_{k=1}^{n/2}f(x_{2k-1}) + 2\sum_{k=1}^{(n/2)-1}f(x_{2k}) + f(x_1)\bigg]\\ h &= \frac{x_1 - x_0}{n}, x_k = x_0 + kh \end{aligned}

# getCoefficient.persudocode

if x == x0 or x == x1: return 1
else if x == x0 + k.odd: return 4
else if x == x0 + k.even: return 2

Algebraic Percision

小毛龙

Lecture 06. Numerical Solutions to Differential Equations

Differential Equations

A differential equation is a relation between the independent variable x, the dependent variable y and its derivatives. Some of these variables might be missing from the equation.

In this course, we only concern about ordinary differential equations (O.D.E.s).

An O.D.E. can be written as:
\begin{aligned} \frac{\text d y}{\text d x} = f(x, y) \end{aligned}

Numerical Solutions

The base of the numerical methods to solve D.E.s is

\begin{aligned} y(x_0+h) \approx y(x_0) + hy'(x_0) \end{aligned}

where y'(x_0) can be approximated by
\begin{aligned} y'(x_0) \approx \frac{y(x_0+h)-y(x_0)}{h} \end{aligned}

Thus a sequence can be written as y_{n+1} \approx y_{n} + hy'_{n}.
In a known initial value problem, term y'_{n} is often written as f(x_n, y_n).

Euler's Method

\begin{aligned} y_{n+1} = y_n + hf(x_n, y_n) \end{aligned}

Backward Euler's Method

Backward Euler's method is a slight modification to the Euler's method, where y'_{n+1} is used instead of y'_{n}.

\begin{aligned} y_{n+1} = y_n + hf(x_{n+1}, y_{n+1}) \end{aligned}

This method is also known as implicit, when using this method, the y_{n+1} appeared on the right needs to be moved to the left of the equation to be solved.

Heun's Method

Heun's method is a modification to the Euler's method, where an average value of slope is used.

\begin{aligned} y_{n+1} &= y_n + h[\frac{1}{2}(f(x_n, y_n) + f(x_{n+1}, y^{*}_{n+1}))]\\ y^{*}_{n+1} &= y_n + hf(x_n, y_n) \end{aligned}

Taylor's Method of Order 2

\begin{aligned} y_{n+1} &\approx y_n + hy'_n +\frac{h^2}{2}y''_n+R_2 \end{aligned}

The R_2 is called the remainder, and can be neglected if h is adequately small.