欧拉方法详解

phenom

高中牛顿力学回顾

有一个具有一定速度在运动的物体:

NIJ9aT.png

当我们需要对其进行模拟时,自然会想起高中的 位移 = 速度 * 时间,即:

$$s = v * t$$

而当该物体具有恒定加速度(恒力)时:

NIJpZV.png

我们可以应用牛顿第二定律,求出物体在 $t$ 时间的速度:

$$v_t = v_0 + \frac Fm * t$$

因此,我们想获得物体在 $t$ 时间的位置,就可以使用以下公式:

$$s = v_0 * t + \frac 12 * \frac Fm * t^2$$

该公式在理想情况下,十分合理且正确,比如模拟一个斜方向向上抛出的小球,重力作用下(无空气摩擦),我们必定能得到一个对称的抛物线:

NIJCIU.png

这都是高中物理课上我们学到的知识,对于简单的模拟已经够用了。

使用数值方法

然而,我们要模拟的现实世界往往比较复杂,几乎不可能有受力恒定的情况出现,考虑以下情况:

NIGxrq.png

你能写出红色盒子在之后的每一个时间 $t$ 的位置方程吗?显然这是无比困难的,因为该盒子受到的外力 $f$ 跟其位置 $p$ 有关,即 $f$ 与 $p$ 存在某种函数关系,因此,我们可以抽象出物体加速度 $a$ 与 $p$ 的关系:

$$a = \frac {f(t, p)} m$$

而又因 $a$ 为 速度 $v$ 的一阶阶导数,$v$ 为 $p$ 的一阶导数,因此我们最终得到:

$$\ddot p = \dot v = \frac {f(t, p)} m$$

这是一个标准的微分方程($p$ 上面有一点表示导数),我们从该方程中直接求出 $p$ 基本上是不可能的,因 $f$ 往往很复杂(包含了各种碰撞和约束) 。但是我们可以使用数值方法数值方法的本质上是用某种方法来近似求出我们需要的解。使用数值方法我们得不到精确解,但是没所谓,在精度要求不是十分严格的模拟下一般都够用了。


而又物体的速度 $v$ 为 $p$ 的一阶导数,我们得到:

$$\dot p = v$$

我们可以将两条微分方程组合起来,使用矩阵表示:

$$\begin{bmatrix}\dot v \\ \dot p\\ \end{bmatrix} = \begin{bmatrix}\frac {f(t, p)} m \\ v\end{bmatrix}$$

用 $\dot X$ 表示左边的矩阵,$F(t, X)$ 表示右边的矩阵,有:

$$\dot X = F(t, X)$$

啰嗦了这么多,就是为了得到这个方程,我们的目标就是要使用数值方法求解出这个方程的 $X$。为什么要求 $X$?因为模拟的本质就是求某个物体在某一时候应该正确出现的位置(和旋转)


求解微分方程的数值方法有许多,如 Runge KuttaMid-point,还有我之前介绍过的 Verlet,但我今天介绍最常用也是最基础的三种方法,分别为:

  • 显式(向前)欧拉方法:Explicit Euler
  • 隐式(向后)欧拉方法:Implicit Euler
  • 半隐式欧拉方法:Semi-implicit Euler

显式(向前)欧拉方法

1. 介绍

观察我们的方程:

$$\dot X = F(t, X)$$

求解 $X$ 不容易,但是我们可以从 $\dot X$ 入手。

回忆高中数学,若有函数 $f(x)$,那么其导数的定义为:

$$\dot f(x) = \frac {df}{dx} = \lim_{h\to 0}\frac {f(x + h) - f(x)}{h}$$

但是计算机模拟是离散的,因此步长 $h$ 取不了无限小,只能是一个有一定大小的数,但是我们可以利用这种思想,对 $\dot f(x)$ 进行近似:

$$\dot f(x) \approx\ \frac {f(x + h) - f(x)}{h}$$

按照这种思路我们可以求出近似的 $\dot X$。

比如,我们要求 $t_1$ 时刻的 $\dot X$,我们可以:

$$\dot X_{1} = F(t_1, X_1) = \frac {X_2 - X_1}{t_2 - t_1}$$

其中 ${t_2 - t_1}$ 即模拟步长 $h$。那么现在我们可以求 $t_2$ 时刻的 $X$,即 $X_2$ 为:

$$X_2 = X_1 + h * F(t_1, X_1)$$


这种对于每一个时刻($t_2$),都使用上一时刻($t_1$)的解进行求解的方法即为显式欧拉方法,又叫向前欧拉方法。下图显示了显式欧拉方法的细节:

NIGvMn.png

绿线为 $\dot X_{1}$ 的真实值,黑线为显式欧拉的得到的近似值,可以看到是有一定误差存在的,同时也能发现,缩短步长 $h$ 能提高精度。

显式欧拉简单高效且容易实现,但是缺点明显。

2. 稳定性分析

显式欧拉方法不是无条件稳定的方法。这句话我们要理解两个词,首先是稳定,第二是无条件。稳定指的是在长时间模拟下,得到的数值不会指数级爆炸;无条件指的是在任何微分方程下都能稳定。

因此显式欧拉方法在某些方程下能稳定,某些方程下不稳定。我们可以用以下方程来测试:

$$\dot y(t) = -\lambda y(t)$$ $$y(0) = \hat y$$

该线性微分方程称为测试方程,其中有初始值 $\hat y$,且 $\lambda > 0$

该方程的解为 $y = \hat y e^{-\lambda t}$,显然在 $t \to +\infty$ 时,$y \to 0$

我们将该测试方程代入到显式欧拉方法的公式里,有:

$$y_{t + 1} = y_t - h\lambda y_t = (1 - h\lambda)y_t $$

运用数学归纳法,我们可以得到 $y_t$ 与 $\hat y$ 的关系为:

$$y_t = {(1 - h\lambda)}^ty_t$$

显然,这是一个指数函数,因此想要在 $t \to +\infty$ 时不发生数值爆炸,我们只能限制:

$$\vert {1 - h\lambda} \vert < 1$$

也就是只有当 $h\lambda > 0$ 时或 $h\lambda < 2$ 时,不等式才能成立。

隐式(向后)欧拉方法

1. 介绍

显式欧拉方法的思路是使用上一时刻的导数($F(t - 1, X_{t - 1})$)去近似当前时刻的结果($X_{t}$)。而隐式欧拉对齐进行了改进:使用当前时刻的导数($F(t, X_{t})$)近似当前时刻的结果($X_t$)。这是隐式欧拉和显式欧拉的唯一区别。用公式表示就是:

$$\dot X_{1} = F(t_1, X_1) = \frac {X_1 - X_0}{t_1 - t_0} $$ $$ \Rightarrow X_1 = X_0 + h * F(t_1, X_1)$$

NIGzq0.png

绿线为 $\dot X_{1}$ 的真实值,黑线为隐式欧拉的得到的近似值,隐式欧拉同样也会产生误差。

也许你会发现一些问题:在 $X_1$ 未求出前,我们怎么能先计算 $F(t_1, X_1)$ 呢?没错这就是隐式欧拉方法最大的缺陷,它计算困难复杂,通常需要求解方程组,但是隐式欧拉最大的优点是数值稳定。

2. 稳定性分析

隐式欧拉方法是无条件稳定的方法。使用上面的测试方程,我们对隐式欧拉方法进行同样的操作:

$$y_{t + 1} = y_t - h\lambda y_{t + 1} $$

整理得:

$$y_{t + 1}(1 + h\lambda) = y\_t$$

对初始值 $\hat y$ 归纳得:

$$y_t = {(\frac {1}{1 + h\lambda})}^t \hat y$$

我们又得到了一个指数函数,同样想要不发生数值爆炸,只有:

$$\vert \frac {1}{1 + h\lambda} \vert < 1$$

也就是需要 $\vert 1 + h\lambda \vert > 1$。步长 $h$ 必定是一个正数,又有 $\lambda > 0$(上面有写),因此不等式恒成立。

半隐式欧拉方法

显式欧拉和隐式欧拉都有各种的缺点,因此后人想出了一种新方法,将两者结合了起来。

现在让我们忘记 $\dot X = F(t, X)$,回到:

$$\dot v = \frac {f(t, p)} m $$ $$ \dot p = v$$

对其使用显式欧拉积分,有:

$$v_{t + 1} = v_t + h * \frac {f(t, p_t)} m $$ $$ p_{t + 1} = p_t + h * v_t$$

既然对两条方程都使用隐式欧拉比较难求解,那么我们就仅对较简单的那条方程使用隐式欧拉:

NIJiiF.png

红框显式了该方法与显式欧拉的唯一不同之处,我们使用显式欧拉求解速度,用隐式欧拉求解位置,这种混合的方法就叫半隐式欧拉方法

半隐式欧拉方法相比于隐式欧拉方法更容易实现,速度更快,相比显式欧拉要更稳定些,但并不是无条件稳定。在大部分情况下半隐式欧拉方法都能稳定,除了某些极端条件(比如非常大的参数)。所以半隐式欧拉在刚体模拟下是最常用,也是首推的。

总结

  • 显式欧拉:简单,快,不稳定
  • 隐式欧拉:复杂,慢,稳定
  • 半隐式欧拉:简单,快,大部分情况稳定(刚体模拟首推)
阅读 6.5k

可视化/计算几何/图形

9 声望
3 粉丝
0 条评论

可视化/计算几何/图形

9 声望
3 粉丝
文章目录
宣传栏