matlab给定状态方程用欧拉法_后向欧拉法——从粒子系统谈起

引言

考虑这样一个常微分问题，我们只知道函数

上的某个点

和该函数的导数

，但我们想要求得函数

的表达式，或者至少要求得函数上的一些离散点，以便我们使用样条曲线或者别的什么方法近似这个函数。

这样的问题其实很常见，以粒子系统来说，我们知道每个粒子的初始位置

和粒子此时的速度和其受到的力，我们想要求得粒子下一个时间点新的位置

。

粒子受到的力包括重力这样的常数力，包括和粒子位置相关的力（比如空间中的力场），包括和粒子速度相关的力（比如空气摩擦力），包括弹力。

一个直观的想法是，我们假设粒子从此刻

时间开始，接下来速度是匀速的。根据小学二年级我们学到的运动学知识，我们能得到经历时间步长

后，粒子的位置为：

这就是赫赫有名的前向欧拉法，你显然知道结果和真实结果存在误差，因为这一个时间步长范围内粒子的速度是会变化的，甚至加速度都会发生变化，我们的假设只有在小学二年级的期末考试试卷上才是合理的。

但是在介绍后向欧拉方法之前，我们先来介绍一下欧拉方法，而这，要从泰勒级数方法说起。

---

从泰勒级数方法说起

不提泰勒级数方法，就直接灌输欧拉方法的行为是完全没有意义的。 ——拜麦瑟夫

我们重新梳理一遍我们的目标，我们需要求得

的函数，或者求得函数上的一些离散点，如下表：

其中,

是

的一个近似值，根据上面的表，可以构造样条函数或者近似函数。而对于粒子系统来说，求这样的表即是我们的目标。

我们的问题表示为：

其中

是一个已知的两个变量的函数，而

为函数上一个已知点，我们需要求出函数

，使得

的某个领域内，有

且

。或者说，给定一个小步长

，我们想要求

。

根据泰勒级数展开可得：

虽然后面还有无穷多项，但我的耐心只有3，所以暂时只写到这里。根据已知条件，我们已经知道了目标函数的一阶导数

，那么我们可以根据此求二阶导数，再求三阶导数，在耐心消失之前尽可能求更多的项，便能近似求出

来。

如果你的耐心和我一样，在计算

中只使用了

以及之前的项，省略号以后的所有项就共同构成了这个方法的

截断误差。

截断误差

因为我们的计算不包含泰勒展开中

以上的项，所以局部截断误差为

，因此，当

的时候局部误差正比于

，如果

，那么

，经历几千次迭代以后就有可能得到完全不正确的解。

许多局部截断误差累积之后引起整体截断误差，如果局部截断误差为

，那么整体截断误差就为

，因为为了到达

位置，需要计算

次。

---

欧拉方法

在开始正式讲解欧拉方法之前，我们回过头来看看泰勒级数方法的优缺点，优点在于足够简单，只需要求导和不停迭代就能求出所需要的值，并且可以有着非常高的精度，前提是你耐心足够。

而缺点是需要反复对

求导，假设这个函数在范围内处处可导，对

关于

求导的话，显然还需要用到链式法则，反复求导势必会耗费大量计算资源。

除非，你只用一阶导，将函数近似为：

不要怀疑，这就是你在本文最开始看到的欧拉方法。你不需要求任何导数，但是由于截断误差非常大你必须选取非常小的

值，导致保证精确度的条件下性能会很差。

后向欧拉法

后向欧拉法的表现形式如下：

$$

如果你的数学感觉和我一样差，你会疑惑为什么这个式子是对的，在此之前的式子就算容易出现大的误差，至少还是一个近似泰勒展开的结果，而这个式子似乎奇奇怪怪。我们先分析一下这个式子的合理性：

已知泰勒展开为：

(1)

我们把

也进行泰勒展开：

把

换到等式左边，再同时乘

:

如果我们把得到的式子带回到 (1)中，就能得到：

按照逐渐失去耐心的老规矩，我们只保留前面的2项，让后面的项默默离开，就能得到后向欧拉法的式子。并且可以证明，后向欧拉法是收敛的，也就是说它的误差不会因为前进多次步长而像前向欧拉法那样发散。

不发散就代表计算结果准确吗？

如果我们已知

的解析表达式的话，似乎没啥缺点，但往往我们无法得到解析表达式（否则为啥还需要数值计算，解析计算它不香吗？），这时候需要解一个方程组。例如导数

不是关于

的函数，而是一个关于

的函数，实际就需要求解：

关于

项，我们还是得用泰勒展开来求，为方便表示，令

：

值得注意的是，这样的表达是一个近似表达，后续还有无穷多项，也就是说，运用欧拉方法求解

的时候，随着迭代，

会不断减小，且步长越大，减小的越快。对于整个系统来说，系统的能量会随着迭代次数而逐渐衰减，这也是常说的隐式欧拉法的damping问题。

这种情况大量出现在粒子系统等物理模拟情况中，每个粒子的受力其实是一个根据位置、速度相关的函数，且难以解析求解。我需要知道下一个时刻的位置和速度才能知道下一个时刻受到的力，但我需要下一个时刻的力才能求出下一个时刻的位置和速度，互相套娃，非常难受。

而随着对力用欧拉方法近似，在整个系统中力就会越来越小，系统能量快速消散。

看到这里，我们已经了解了解常微分方程的欧拉方法，了解了截断误差，知道前向欧拉法为什么算着算着结果就不受控制，知道了后向欧拉法相比起前向欧拉法来说damp会更快，也知道如果需要求解粒子系统，我们需要求解一个方程组。

那么，让我们开始，进入正题吧！

---

求解粒子系统

首先明确我们需要求解怎样的方程组。

我们知道力是关于位置，速度的函数：

，我们需要计算的是

时刻的

和

，根据后向欧拉法，我们能得到的方程组为：

中学的物理知识告诉我们，速度是位置的一阶导，加速度是速度的一阶导，把(1)式稍稍变形所以

为了简化表示，令

将

$泰勒展开：

这里的

都是偏导。

把展开后的式子中

通过(3)式替换成

,统一带回(2)，则得到一个只关于速度的式子：

我们知道

关于

的偏导，知道当前的速度和位置，上述方程只有

为未知数，当求得 后，我们自然求得下一个状态的速度，也即求出了下一个状态位置的导数，即求算出最终的位置和速度。