matlab给定状态⽅程⽤欧拉法_后向欧拉法——从粒⼦系统谈
起
引⾔
考虑这样⼀个常微分问题，我们只知道函数
上的某个点
和该函数的导数
，但我们想要求得函数
的表达式，或者⾄少要求得函数上的⼀些离散点，以便我们使⽤样条曲线或者别的什么⽅法近似这个函数。
这样的问题其实很常见，以粒⼦系统来说，我们知道每个粒⼦的初始位置
和粒⼦此时的速度和其受到的⼒，我们想要求得粒⼦下⼀个时间点新的位置
。
粒⼦受到的⼒包括重⼒这样的常数⼒，包括和粒⼦位置相关的⼒（⽐如空间中的⼒场），包括和粒⼦速度相关的⼒（⽐如空⽓摩擦⼒），包括弹⼒。
⼀个直观的想法是，我们假设粒⼦从此刻
时间开始，接下来速度是匀速的。根据⼩学⼆年级我们学到的运动学知识，我们能得到经历时间步长
后，粒⼦的位置为：
这就是赫赫有名的前向欧拉法，你显然知道结果和真实结果存在误差，因为这⼀个时间步长范围内粒⼦的速度是会变化的，甚⾄加速度都会发⽣变化，我们的假设只有在⼩学⼆年级的期末考试试卷上才是合理的。
但是在介绍后向欧拉⽅法之前，我们先来介绍⼀下欧拉⽅法，⽽这，要从泰勒级数⽅法说起。
从泰勒级数⽅法说起
不提泰勒级数⽅法，就直接灌输欧拉⽅法的⾏为是完全没有意义的。 ——拜麦瑟夫
我们重新梳理⼀遍我们的⽬标，我们需要求得
的函数，或者求得函数上的⼀些离散点，如下表：
其中,
是
的⼀个近似值，根据上⾯的表，可以构造样条函数或者近似函数。⽽对于粒⼦系统来说，求这样的表即是我们的⽬标。
我们的问题表⽰为：
其中
是⼀个已知的两个变量的函数，⽽
为函数上⼀个已知点，我们需要求出函数
，使得
的某个领域内，有
且
。或者说，给定⼀个⼩步长
，我们想要求
。
根据泰勒级数展开可得：
虽然后⾯还有⽆穷多项，但我的耐⼼只有3，所以暂时只写到这⾥。根据已知条件，我们已经知道了⽬标函数的⼀阶导数
，那么我们可以根据此求⼆阶导数，再求三阶导数，在耐⼼消失之前尽可能求更多的项，便能近似求出
来。
如果你的耐⼼和我⼀样，在计算
中只使⽤了
以及之前的项，省略号以后的所有项就共同构成了这个⽅法的
截断误差。
截断误差
截断误差
因为我们的计算不包含泰勒展开中
以上的项，所以局部截断误差为
，因此，当
的时候局部误差正⽐于
，如果
，那么
，经历⼏千次迭代以后就有可能得到完全不正确的解。
许多局部截断误差累积之后引起整体截断误差，如果局部截断误差为
，那么整体截断误差就为
，因为为了到达
位置，需要计算
次。
欧拉⽅法
在开始正式讲解欧拉⽅法之前，我们回过头来看看泰勒级数⽅法的优缺点，优点在于⾜够简单，只需要求导和不停迭代就能求出所需要的值，并且可以有着⾮常⾼的精度，前提是你耐⼼⾜够。
⽽缺点是需要反复对
求导，假设这个函数在范围内处处可导，对
关于
求导的话，显然还需要⽤到链式法则，反复求导势必会耗费⼤量计算资源。
除⾮，你只⽤⼀阶导，将函数近似为：
不要怀疑，这就是你在本⽂最开始看到的欧拉⽅法。你不需要求任何导数，但是由于截断误差⾮常⼤你必须选取⾮常⼩的
值，导致保证精确度的条件下性能会很差。
后向欧拉法
后向欧拉法的表现形式如下：
$$
如果你的数学感觉和我⼀样差，你会疑惑为什么这个式⼦是对的，在此之前的式⼦就算容易出现⼤的误差，⾄少还是⼀个近似泰勒展开的结果，⽽这个式⼦似乎奇奇怪怪。我们先分析⼀下这个式⼦的合理性：
已知泰勒展开为：
(1)
我们把
也进⾏泰勒展开：
把
换到等式左边，再同时乘
:
如果我们把得到的式⼦带回到 (1)中，就能得到：
按照逐渐失去耐⼼的⽼规矩，我们只保留前⾯的2项，让后⾯的项默默离开，就能得到后向欧拉法的式⼦。并且可以证明，后向欧拉法是收敛的，也就是说它的误差不会因为前进多次步长⽽像前向欧拉法那样发散。
不发散就代表计算结果准确吗？
如果我们已知
的解析表达式的话，似乎没啥缺点，但往往我们⽆法得到解析表达式（否则为啥还需要数值计算，解析计算它不⾹吗？），这时候需要解
matlab求导⼀个⽅程组。例如导数
不是关于
的函数，⽽是⼀个关于
的函数，实际就需要求解：
关于
项，我们还是得⽤泰勒展开来求，为⽅便表⽰，令
：
值得注意的是，这样的表达是⼀个近似表达，后续还有⽆穷多项，也就是说，运⽤欧拉⽅法求解
的时候，随着迭代，
会不断减⼩，且步长越⼤，减⼩的越快。对于整个系统来说，系统的能量会随着迭代次数⽽逐渐衰减，这也是常说的隐式欧拉法的
damping问题。
这种情况⼤量出现在粒⼦系统等物理模拟情况中，每个粒⼦的受⼒其实是⼀个根据位置、速度相关的函数，且难以解析求解。我需要知道下⼀个时刻的位置和速度才能知道下⼀个时刻受到的⼒，但我需
要下⼀个时刻的⼒才能求出下⼀个时刻的位置和速度，互相套娃，⾮常难受。
⽽随着对⼒⽤欧拉⽅法近似，在整个系统中⼒就会越来越⼩，系统能量快速消散。
看到这⾥，我们已经了解了解常微分⽅程的欧拉⽅法，了解了截断误差，知道前向欧拉法为什么算着算着结果就不受控制，知道了后向欧拉法相⽐起前向欧拉法来说damp会更快，也知道如果需要求解粒⼦系统，我们需要求解⼀个⽅程组。
那么，让我们开始，进⼊正题吧！
求解粒⼦系统
⾸先明确我们需要求解怎样的⽅程组。
我们知道⼒是关于位置，速度的函数：
，我们需要计算的是
时刻的
和
，根据后向欧拉法，我们能得到的⽅程组为：
中学的物理知识告诉我们，速度是位置的⼀阶导，加速度是速度的⼀阶导，把(1)式稍稍变形所以
为了简化表⽰，令
将
$泰勒展开：
这⾥的
都是偏导。
把展开后的式⼦中
通过(3)式替换成
，统⼀带回(2)，则得到⼀个只关于速度的式⼦：
我们知道
关于
的偏导，知道当前的速度和位置，上述⽅程只有
为未知数，当求得 后，我们⾃然求得下⼀个状态的速度，也即求出了下⼀个状态位置的导数，即求算出最终的位置和速度。