欧拉法是一种求解微分方程的数值方法。微分方程的目标是求出方程的具体表达式，但数值方法则是根据方程求出其每一个点，而不是找出表达式，在工程上具有很大的实用性。

欧拉法求解微分方程组的原理

设有微分方程：
$$\frac{dx(t)}{dt}=f(x)$$
$x(t_0)=x_0$已知。

我们对上述方程进行积分，积分区域为$[t_0,t_1],\Delta t=t_1-t_0$，可得：
$$x(t_1)=x(t_0)+{\int }_{t_0}^{t_1}f(t)dt$$
其中$x(t_1)$就是我们要求解的东西了。于是，就差积分怎么解决了。

对于欧拉法来说，积分是靠梯形面积来近似，如下所示：
$${\int }_{t_0}^{t_1}f(t)dt=f(x(t_0))(t_1-t_0)=f(x_0)\Delta t$$
带入上述方程即可得：
$$x(t_1)=x(t_0)+f(x_0)\Delta t$$
按照上式进行递推，即可得：
$$x_{k+1}=x_k+f(x_k)\Delta t$$
其中 $x_0$ 已知。问题就搞定了。

*欧拉法原理图解

因为
$$f(x_k)=\frac{dx}{dt}{|}_{t=t_k}$$
所以递推式可以写成：
$$x_{k+1}=x_k+\frac{dx}{dt}{|}_{t=t_k}\Delta t$$
而$\frac{dx}{dt}$是所求方程$x(t)$的斜率，因此实际上是将下图的阴影部分的面积来替代积分表达式：

欧拉法数值稳定性讨论

方便讨论，我们取一种简单情况，设有微分方程：
$$\frac{dx(t)}{dt}=ax$$
什么是数值稳定性问题?

假设 k 步计算中，$x_k$与实际值存在误差${\rho }_k$，则第 k+1 步计算中，误差是否会增大。若是，则不稳定啦。因为越往后计算，这个误差肯定是越爆炸的。

那么欧拉法数值稳定性如何呢？首先假设$x_k$与实际值存在误差${\rho }_k$，根据递推公式：
$$x_{k+1}=x_k+f(x_k)\Delta t$$
设$f(x_k)=a{x}_k$而：
$${\hat{x}}_k=x_k+{\rho }_k$$
故第 k+1 步时时：
$${\hat{x}}_k=(1+a\Delta t)(x_k+{\rho }_k)=x_{k+1}+(1+a\Delta t){\rho }_k$$
于是 k+1 步的误差为：
$${\rho }_{k+1}=(1+a\Delta t){\rho }_k$$
若要求稳定，则必须有$|1+a\Delta t|\le 1$
于是可得欧拉法的数值稳定区域如下：

欧拉法误差分析——局部截断误差

微分方程$\frac{dx}{dt}=f(x)$，用欧拉法时，我们的迭代方程是：
$x_k=x_{k-1}+f(x_{k-1})\Delta t=x_{k-1}+\frac{dx}{dt}{|}_{t_{k-1}}\Delta t$。

若对 $x_k$ 进行泰勒展开，可得：
$$x(t_k)=x(t_{k-1}+\Delta t)=x_{k-1}+x^{\prime }(t_{k-1}\Delta t+x^{\prime \prime }(t_{k-1})\Delta t^2+\cdots$$
因此，欧拉法的截断误差为：
$$x^{\prime \prime }(t_{k-1})\Delta t^2+\cdots$$
由于有$\Delta t^2$，所以我称欧拉法的精度为 2-1 = 1 阶。

改进欧拉法原理

欧拉法是用梯形的面积来近似积分的，而且梯形的上底是用 KaTeX parse error: Expected 'EOF', got '}' at position 7: x_[i-1}̲时，曲线$x(t)$的斜率来近似计算的，所以不太好，可以看到的是，欧拉法的精度只有一阶，其阶段误差与 $\Delta t^2$有关。这个误差之所以那么大，是因为我们在求积分（面积）的时候，用了以$x_{k-1}$为上底的梯形去近似他。因此，我们需要对这个梯形的上底进行一定的修正。

首先，根据递推公式求出：
$$x_k^{(0)}=x_{k-1}+f(x_{k-1})\Delta t$$
之后再求出修正后的 $x_k$：
$$x_k=x_{k-1}+\frac{1}{2}(f(x_{k-1}+x_k^{(0)})\Delta t$$
其原理图如下：

改进欧拉法的稳定性分析

设微分方程为 $\frac{dx}{dt}=\lambda x$，于是KaTeX parse error: Undefined control sequence: \ambda at position 16: x_{k}^{(0)]=(1+\̲a̲m̲b̲d̲a̲\Delta t)x_{k-1…，$f(x_k^{(0)}=\lambda (1+\lambda \Delta t)x_{k-1}$，于是$x_k=(1+\lambda \Delta t+\frac{1}{2}(\lambda \Delta t{)}^2)x_{k-1}$

所以误差不增值的条件为：$|1+\lambda \Delta t+\frac{1}{2}(\lambda \Delta t{)}^2|\le 1$，表现在图上，即为：

对比欧拉法，显然稳定区域大了一点。

改进欧拉法的误差分析

因为 $f(x_k^0)=\frac{dx}{dt}{|}_{x_{k-1}}+\frac{d^{2}x}{d{t}^2}{|}_{x_{k-1}}\Delta t^2+O(\Delta t^2)$

于是$$x_k=x_{k-1}+\frac{dx}{dt}{|}_{x_{k-1}}\Delta t+\frac{d^{2}x}{d{t}^2}{|}_{x_{k-1}}\Delta t^2$$
对比泰勒展开，可知其误差为$O(\Delta t^3)$，因此精确度为 3-1 = 2阶！

隐式梯形法

隐式梯形法也是利用 $x_k$ 来改善上底，只是其递推公式为：
$$x_k=x_{k-1}+\frac{\Delta t}{2}(f(x_{k-1})+f(x_k))$$
因此，上底的修正，包含在上述方程中。因此，不同于欧拉法，每一次迭代，隐式梯形法都要求解一次代数方程。

隐式梯形法的稳定性和误差

首先是稳定性，同样对于微分方程 $\frac{dx}{dt}=\lambda x$，可得其误差的表达式为：

误差不增值条件:

整理可得：
因此其稳定区域为：复平面的左半部分，可见比改进欧拉发、欧拉法要大得多！！

至于误差，同样可以证明的是隐式梯形法的精确度为2阶。

案例

可见求解区域为 [0，1]，我们设置求解的步数为 100，也即 $\Delta t=0.01s$，代码如下：

import numpy as np
import matplotlib.pyplot as plt
def f(x):
    return -20*x

def uler(n,x0,t0,tn,f):
    x = [x0]
    t = [t0]
    for i in range(n):
        dt = (tn-t0)/n
        tk = t[i]+dt
        t.append(tk)
        xk = x[i]+f(x[i])*dt
        x.append(xk)
    return x,t

def uler2(n,x0,t0,tn,f):
    x = [x0]
    t = [t0]
    for i in range(n):
        dt = (tn-t0)/n
        tk = t[i]+dt
        t.append(tk)
        xkk = x[i]+f(x[i])*dt
        xk = x[i]+ 1/2*(f(x[i])+f(xkk))*dt
        x.append(xk)
    return x,t


import sympy
x = sympy.symbols('x')

def uler3(n,x0,t0,tn,f):
    x = [x0]
    t = [t0]
    for i in range(n):
        dt = (tn-t0)/n
        tk = t[i]+dt
        t.append(tk)
        xk = (x[i]+dt*f(x[i])/2)/(1+10*dt)
        x.append(xk)
    return x,t
        
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False

font1 = {
    'family' : 'SimHei',
'weight' : 'normal',
'size'   : 15,
}

def myplot(n,x0,t0,tn,f):
    x1,t = uler(n,x0,t0,tn,f)
    x2,t = uler2(n,x0,t0,tn,f)
    plt.figure(figsize=(12,4))
    plt.subplot(1,3,1)
    plt.plot(t,x1,linewidth=3,color='r',label='欧拉法求解')    
    plt.xlabel('t',fontsize=24)
    plt.ylabel('x',fontsize=24)
    plt.legend(prop=font1)
    plt.grid()
    
    plt.subplot(1,3,2)
    plt.plot(t,x2,linewidth=3,color='b',label='改进欧拉法求解')
    plt.xlabel('t',fontsize=24)
    plt.ylabel('x',fontsize=24)
    plt.legend(prop=font1)
    plt.grid()
    
    x3,t = uler3(n,x0,t0,tn,f)
    plt.subplot(1,3,3)
    plt.plot(t,x3,linewidth=3,color='g',label='隐式梯形法求解')
    plt.xlabel('t',fontsize=24)
    plt.ylabel('x',fontsize=24)
    plt.legend(prop=font1)
    plt.grid()
    
if __name__ == '__main__':
    myplot(100,1,0,1,f)

求解结果：

可以修改步长为 15，可得：

为什么步数为 15，结果就变形了呢？这是因为算法的稳定性、误差的原因，这里不再分析。