Common Numerical Integration and Linear System Solving Methods

Numerical integration

For a certain function $f(x)$, we need to calculate the future state of $t+h$ from the state at time t. Calculate $f(t+h)$ from $f(t)$. To this end, we will expand $f(t)$ by Taylor expansion $f(x+\mathrm{\Delta }x)=f(x)+f^{\prime }(x)\mathrm{\Delta }x+\frac{1}{2}{f}^{″}(x)\mathrm{\Delta }{x}^2+\frac{1}{6}{f}^{‴}(x)\mathrm{\Delta }{x}^3+\cdots$, so the following methods are all approximate truncations.We can classify the integration methods by taking the highest order of the Taylor expansion. If Euler integral only uses the first derivative, it is a first-order method. And its error is $O(n^2)$, and after halving the step size $h$, the error becomes a quarter of the original error.

Euler integral

The classical expression of approaching curves with a line

Explicit Euler

The simplest integration method, however, may have numerical expansion

$\overrightarrow{v}(t_1)=\overrightarrow{v}(t_0)+M^{-1}\overrightarrow{F}(t_0)\mathrm{\Delta }t$

$x(t_1)=x(t_0)+\overrightarrow{v}(t_0)\mathrm{\Delta }t$

Implicit Euler

Updating the present with future states requires solving linear systems with numerical dissipation

$\overrightarrow{v}(t_1)=\overrightarrow{v}(t_0)+M^{-1}\overrightarrow{F}(t_1)\mathrm{\Delta }t$

$x(t_1)=x(t_0)+\overrightarrow{v}(t_1)\mathrm{\Delta }t$

Semi_Implicit Euler

The combination of explicit and implicit Euler for numerical stability

$\overrightarrow{v}(t_1)=\overrightarrow{v}(t_0)+M^{-1}\overrightarrow{F}(t_0)\mathrm{\Delta }t$

$x(t_1)=x(t_0)+\overrightarrow{v}(t_1)\mathrm{\Delta }t$

Leap Frog

A variant of semi implicit Euler, which separates velocity and displacement by half an hour and updates them using the midpoint method, resulting in smaller errors

$\overrightarrow{v}(t_{0.5})=\overrightarrow{v}(t_{-0.5})+M^{-1}\overrightarrow{F}(t_0)\mathrm{\Delta }t$

$x(t_1)=x(t_0)+\overrightarrow{v}(t_{0.5})\mathrm{\Delta }t$

Velocity verlet

Introducing velocity to simplify the Welley integral and changing the midpoint method to the trapezoidal method compared to the frog jumping method, so that velocity and position can be synchronized in time

$x(t_1)=x(t_0)+\overrightarrow{v}(t_0)\mathrm{\Delta }t+\frac{1}{2}{M}^{-1}\overrightarrow{F}(t_0)\mathrm{\Delta }{t}^2$

$\overrightarrow{v}(t_1)=\overrightarrow{v}(t_0)+\frac{1}{2}{M}^{-1}(\overrightarrow{F}(t_0)+\overrightarrow{F}(t_1))\mathrm{\Delta }t$

Runge Kutta

The Euler method is actually a first-order Runge Kutta, while the higher-order Runge Kutta uses an approximation of the derivative at the midpoint to estimate the subsequent derivative using the previous derivative

RK2

$\overrightarrow{v}(t_1)=\overrightarrow{v}(t_0)+M^{-1}\overrightarrow{F}(t_0)\mathrm{\Delta }t$

$x(t_{0.5})=x(t_0)+\frac{1}{2}\overrightarrow{v}(t_0)\mathrm{\Delta }t$

$\overrightarrow{v}(t_{0.5})=\overrightarrow{v}(t_0)+\frac{1}{2}{M}^{-1}\overrightarrow{F}(t_{0.5})\mathrm{\Delta }t$

$x(t_1)=x(t_0)+\overrightarrow{v}(t_{0.5})\mathrm{\Delta }t$

Linear system

The classical solution to the problem of $Ax=b$ is too costly to directly solve the inverse, so each individual needs to optimize based on their own abilities

Matrix decomposition

Decompose a complex matrix into the product of several simple matrices to facilitate inversion. Taking the eigen library as an example, the most commonly used ones are LU decomposition (which decomposes into the product of lower triangular matrix $L$ and upper triangular matrix $U$), LLT decomposition (which requires positive definite symmetry of $A$, decomposes into the product of lower triangular matrix $L$ and its transposed $L^T$, with a speed twice that of LU decomposition), and the specific method is Gaussian elimination

Linear iteration

Stop after a certain number of iterations or reaching a certain accuracy

Jacobi iteration

For $Ax=b$, splitting $A$ into $D-L-R$
$$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]b=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]D=\left[\begin{array}{ccc}{a}_{11}& 0& 0\\ 0& a_{22}& 0\\ 0& 0& a_{33}\end{array}\right]$$

$$L=\left[\begin{array}{ccc}0& 0& 0\\ a_{21}& 0& 0\\ a_{31}& a_{32}& 0\end{array}\right]R=\left[\begin{array}{ccc}0& a_{12}& a_{13}\\ 0& 0& a_{23}\\ 0& 0& 0\end{array}\right]f=\left[\begin{array}{c}\frac{b_1}{a_{11}}\\ \frac{b_2}{a_{22}}\\ \frac{b_3}{a_{33}}\end{array}\right]$$
$$B=I-D^{-1}A，x_{k+1}=B{x}_k+f$$

Gauss-Seidel iteration

In Jacobi iteration, the result will only be applied to the next round after each round of calculation is completed. If it is updated in a timely manner after the current round of calculation, the efficiency is intuitively better than Jacobian iteration, but Jacobian iteration is naturally suitable for parallelism
$$Jacobi\{\begin{array}{l}\begin{array}{c}{x}_1^{k+1}=\frac{1}{a_{11}}(-a_{12}{x}_2^k-a_{13}{x}_3^k+b_1)\\ x_2^{k+1}=\frac{1}{a_{22}}(-a_{21}{x}_1^k-a_{23}{x}_3^k+b_2)\\ x_3^{k+1}=\frac{1}{a_{33}}(-a_{31}{x}_1^k-a_{32}{x}_2^k+b_3)\end{array}\end{array}$$

$$Gauss-Seidel\{\begin{array}{l}\begin{array}{c}{x}_1^{k+1}=\frac{1}{a_{11}}(-a_{12}{x}_2^k-a_{13}{x}_3^k+b_1)\\ x_2^{k+1}=\frac{1}{a_{22}}(-a_{21}{x}_1^{k+1}-a_{23}{x}_3^k+b_2)\\ x_3^{k+1}=\frac{1}{a_{33}}(-a_{31}{x}_1^{k+1}-a_{32}{x}_2^{k+1}+b_3)\end{array}\end{array}$$