数值积分

对于某一函数$f(x)$，我们需要由 $t$ 时刻的状态求得未来某时间 $t+h$ 的状态。即由 $f(t)$求得 $f(t+h)$ 。为此，我们将$f(t)$ 泰勒展开 $f(x+\Delta x)=f(x)+f’(x)\Delta x+\frac{1}{2}f’’(x) \Delta x^2+\frac{1}{6}f’’’(x) \Delta x^3+\cdots$ 实际中不可能取到无穷项，故以下方法均为其的近似截断。

我们可以通过所取泰勒展开中最高阶数对积分方式进行分类。如欧拉积分只用到了一阶导数，所以是一阶方法。且其误差为 $O(n^2)$ ，步长 $h$ 减半后，误差变为原误差的四分之一。

欧拉积分

经典以折线逼近曲线的趋近式表达

显式欧拉(Explicit Euler)

最简单的积分方式，然而会有数值膨胀

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

$x(t_1)=x(t_0)+\vec v(t_0)\Delta t$

隐式欧拉(Implicit Euler)

以未来的状态更新现在，需要解线性系统，且有数值耗散

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

$ x(t_1)=x(t_0)+\vec v(t_1)\Delta t $

半隐式欧拉(Semi_Implicit Euler)

显隐式欧拉的结合，数值稳定

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

$x(t_1)=x(t_0)+\vec v(t_1)\Delta t$

跳蛙法(Leap Frog)

半隐式欧拉的变体，将速度和位移错开半个时步，用中点法更新，误差更小

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

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

速度韦尔莱（Velocity verlet）

引入速度简化韦尔莱积分，与跳蛙法相比改中点法为梯形法，使速度和位置能在时间上同步

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

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

龙格库塔(Runge Kutta)

欧拉法其实就是一阶龙格库塔，高阶龙格库塔则是使用中点处导数的近似值，用之前的导数估算后面的导数

二阶龙格库塔(RK2)

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

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

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

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

线性系统

经典解 $Ax=b$问题，直接求逆过于费，所以需要各凭本事进行优化

矩阵分解

将一个复杂的矩阵拆解为数个简单的矩阵的乘积，从而方便求逆。这里以eigen库为例，最常用的有LU分解（分解为下三角阵 $L$和上三角阵 $U$ 的积），LLT分解（要求 $A$ 正定对称，分解为下三角阵 $L$ 和其转置 $L^T$ 的积，速度两倍于LU分解），具体方法为高斯消元

线性迭代

迭代一定次数或到达一定精度后即可停止

雅可比迭代(Jacobi)

对$Ax=b$，裂解$A$为$D-L-R$
$$
A= \left[ \begin{matrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{matrix} \right] \ b= \left[ \begin{matrix} b_1 \\ b_2 \\ b_3 \end{matrix} \right] \ D= \left[ \begin{matrix} a_{11} & 0 & 0 \\ 0 & a_{22} & 0 \\ 0 & 0 & a_{33} \end{matrix} \right]
$$

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

高斯赛德尔迭代(Gauss-Seidel)

雅可比迭代中每计算完一轮，才会将结果作用于下一轮。若是在当轮计算后就及时更新，直观上效率优于雅可比迭代，但雅可比迭代天然适于并行
$$
雅可比 \begin{cases} \begin{matrix} 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{matrix} \end{cases}
$$

$$
高斯赛德尔\begin{cases} \begin{matrix} 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{matrix} \end{cases}
$$