Summary of numerical analysis algorithm

A summary of the algorithms for numerical analysis, with a brief description of the various methods in Python. Prepare for the exam.

The code given in this paper is mainly for closed-book exam back algorithm. I have to memorize mathematical formulas, just like writing for LaTeX, which is a much scarier thing than writing code. So, some of the main algorithms in the program to write out, easy to remember.

(One part is written when reviewing before the exam, through the sampling inspection of the examination room, more reliable. But at that time, I wrote a little bit more, and then I added a little bit. At this time, I have finished my grades and forgotten everything I learned, so it may not be right. Anyway, don’t expect too much.

If you need more complete code implementation and description of various algorithms, please go to:

  • Github.com/cdfmlr/Nume…
  • The or gitee.com/cdfmlr/Nume…

All right, talk is cheap, let me show you the code!

Take roots of nonlinear equations

dichotomy

x = (a + b) / 2
while True:
    if f(x) * f(a) < 0:
        b = x
    else:
		a = x
	x = (a + b) / 2
Copy the code

Fixed point iteration

x = x_0
while True:
    x = phi(x)
Copy the code

Convergence: Diff (phi, x) exists and is continuous && abs(diff(phi, x)(x_0)) < 1.

Or: If x∗x^*x∗ is the root of f(x)=0f(x)=0f(x)=0, then x∗=… ∗ (x) x ^ * =… X (x ^ *) ∗ =… (x∗), so x∗x^*x∗ is the fixed point of x = phi(x), for x∗x^*x∗ field x: ABS (diff(phi, x)) < 1

Newton iteration

x = x_0
while True:
    x -= f(x) / df(x)
Copy the code

The interpolation

Though laser interpolation

def lagrange_interpolate(points) :
	L = 0  # Interpolation polynomial
    for i, (xi, yi) in enumerate(points):
        li = 1
        for j, (xj, yj) in enumerate(points):
            if j == i:
                continue
            li *= (x - xj) / (xi - xj)
        L += yi * li
    return L
Copy the code

This program is not very easy to understand, or look at the formula:


L ( x ) = ∑ j = 0 k y j ℓ j ( x ) L(x)=\sum _{j=0}^{k}y_{j}\ell _{j}(x)

The P.S. formula is a direct copy of Wikipedia, but KaTeX can’t render it, so it’s a picture. Here is the source code:

$$
\ell _{j}(x)=\prod _{{i=0,\,i\neq j}}^{{k}}{\frac {x-x_{i}}{x_{j}-x_{i}}}={\frac {(x-x_{0})}{(x_{j}-x_{0})}}\cdots {\frac {(x-x_{{j-1}})}{(x_{j}-x_{{j-1}})}}{\frac {(x-x_{{j+1}})}{(x_{j}-x_{{j+1}})}}\cdots {\frac {(x-x_{{k}})}{(x_{j}-x_{{k}})}}
$$
Copy the code

Newton interpolation

Difference quotient:

def dq(f, xs) :
    if len(xs) == 1:  Order # 0
        return f(xs[0])
    # n order
    return (dq(f, xs[1:]) - dq(f, xs[:-1])) / (xs[-1] - xs[0]) 
Copy the code

Hand count table:

Newton interpolation:

def newton_interpolate(points) :
    N = 0
    xs = []
    for (x, y) in points:
        xs.append(x)
        N += dq(f, xs) * prod( [('x' - xi) for xi in xs[:-1]])return N
Copy the code

Least squares fitting

x = [0.1.2.3].T
y = [1.2.4.8].T

X = [1, x, x**2]
# (X.T @ X) @ theta = X.T @ y
theta = pinv(X.T @ X) @ X.T @ y
Copy the code

@ is matrix multiplication (this is the standard Python operator: docs.python.org/3/library/o…

Numerical integration

Trapezoidal quadrature formula

df = (b - a) * (f(a) + f(b)) / 2
Copy the code

Simpson’s quadrature formula

df = (b - a) * (f(a) + 4 * f((a + b) / 2) + f(b)) / 6
Copy the code

Newton-Cotes

Ck(n)C_k^{(n)}Ck(n):

def costes_coefficient(n, k) :
    ckn = ((-1) ** (n - k)) / n * factorial(k) * factorial(n - k)

    h = 1
    t = Symbol('t')
    for j in range(n+1) :ifj ! = k: h *= (t - j) ckn *= integrate(h, (t,0, n))

    return ckn
Copy the code

Make a table to calculate by hand:

Newton-cortez quadrature formula:

def newton_cotes_integral(f, a, b, n) :
    step = (b - a) / n
    xs = [a + i * step for i in range(n+1)]
    return (b - a) * sum([costes_coefficient(n, k) * f(xs[k]) for k in range(0, n+1)])
Copy the code

Initial value problems for ordinary differential equations

Question:


{ y ‘ ( x ) = f ( x . y ) y ( a ) = y 0 ( a Or less x Or less b ) \begin{cases} y'(x)=f(x,y)\\ y(a)=y_0\\ \end{cases} \qquad (a\le x \le b)

Improved Euler

def improved_euler(f, a, b, h, y0) :    
    x = a
    y = y0

    while x <= b:
        yield (x, y)

        y_next_g = y + h * f(x, y)  # forecast
        y_next = y + h * ( f(x, y) + f(x+h, y_next_g) ) / 2  # correction
        
        x = x + h
        y = y_next
Copy the code

RK4

def runge_kutta(f, a, b, h, y0) :
    x = a
    y = y0

    while x <= b:
        yield (x, y)

        k1 = f(x, y)
        k2 = f(x + h / 2, y + h * k1 / 2)
        k3 = f(x + h / 2, y + h * k2 / 2)
        k4 = f(x + h, y + h * k3)
        
        x = x + h
        y = y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
Copy the code

Direct solution of linear equations

Gaussian elimination

A = np.c_[a, b]  # Augmented matrix

n = A.shape[0]
x = np.zeros(n)  # solution

# elimination
for k in range(n-1) :Select pivot entries for the # column. If you do sequential elimination, you don't have to do the next two rows
    i_max = k + argmax(abs(A[k:n, k]))
    A[[i_max, k]] = A[[k, i_max]]
    
    # if A[k][k] == 0
    for i in range(k+1, n):
        m = A[i][k] / A[k][k];
        A[i][k] = 0;
        for j in range(k+1, n+1):
            A[i][j] -= A[k][j] * m

# Back substitution: solve trigonometric equations
x[n-1] = A[n-1][n] / A[n-1][n-1]
for k in range(n-2, -1, -1) :# from n-2 (included) to 0 (included)
    for j in range(k+1, n):
        A[k][n] -= A[k][j] * x[j]
    x[k] = A[k][n] / A[k][k]
Copy the code

LU decomposition

Coefficient matrix A:

n = a.shape[0]
p = [0.1. , n-1]  # Record exchange

for k in range(n-1) :ifI_max = k + argmax(argmax)abs(a[k:n, k]))
        ifi_max ! = k: a[[i_max, k]] = a[[k, i_max]]# swap rows
            p[[i_max, k]] = p[[k, i_max]]  # record

    asserta[k][k] ! =0."Error: pivot is zero"

    for i in range(k+1, n):
        a[i][k] /= a[k][k]  # L @ Strictly lower triangle
        for j in range(k+1, n):
            a[i][j] -= a[i][k] * a[k][j]  # U @ upper triangle
Copy the code

Make the same row swap (p record) for the right end constant B:

b = [b[v] for v in p]
Copy the code

Then we can solve the equation:

Equation ("L * y = b"=> y => y"U * x = y") => x
Copy the code

Add: To solve trigonometric equations:

x = np.zeros(n, dtype=np.float)
x[n-1] = y[n-1] / u[n-1][n-1]
for i in range(n-2, -1, -1) :# from n-2 (included) to 0 (included)
    yi = y[i]
    for j in range(i+1, n):
        yi -= x[j] * u[i][j]
    x[i] = yi / u[i][i]
Copy the code

Iterative solution of linear equations

Jacobi iteration

# A = D - L - U
D = diag(diag(A))
L = - tril(A, -1)
U = - triu(A, 1)

B = inv(D) @ (L + U)
f = inv(D) @ 

while True:
    x_prev = x
    x = B @ x + f
Copy the code

Emmmm, that’s not realistic. I vaguely remember that it would have been easier to just write the equation (I’ve been taking the exam so long I’ve forgotten)…

Gauss – Seidel iterative

B = inv(D - L) @ U
f = inv(D - L) @ b
Copy the code

Eigenvalue method

Is a power law

A = array(shape=(n, n)) # A is an n by n matrix that requires eigenvalues

m = m0 = 1  # Principal eigenvalues
u = u0 = [1.1. .1]  The eigenvector corresponding to #

while True:
	m_prev = m

	v = dot(A, u)
	m = v[argmax(abs(v))]
	u = v / m
Copy the code

The inverse power method

The positive power method is changed:

v = solve(A, u)
Copy the code

So it’s going to be 1 over m and u.

EOF

By CDFMLR 2020.12.31