数值解常微分方程(入门)

                     

贡献者: addis

  1前面几节,我们把求解线性方程组的基本数值方法做了详细的介绍和分析。从这一节开始,我们来尝试使用这些解法,处理更复杂也更贴近实际应用的问题。需要复习小伙伴,可以先去看看前面的内容。

   这一节我们来讨论一下常微分方程(Ordinary Differential Equation(ODE))。注意,我的整个专栏都是在讨论关于科学计算和数值分析的内容。对于常微分方程的分析特性,解析特解和通解等特性,请参考数学分析课程等。

举个例子

   比如,$y'(t)= \cos\left(yt\right) $,这实际上是一个非线性的常微分方程,而且 $y$ 还隐含的是 $t$ 的函数。

   再比如,$y'(t)=ay(t)+b$,如果 $a$ 和 $b$ 都是常数的话,这就是一个线性常微分方程。

   再再比如,$y'(t)=at+b$,求解这个方程可以轻松的将右边对 $t$ 做积分。

   当然,我们也可以写出常微分方程的标准形式

\begin{equation} y'(t)=f(y(t),t),t _0 < t < T,\quad y(t_0)=\hat{y}~. \end{equation}
这里 $y(t_0)=\hat{y}$ 是方程的初始条件,因此这样的问题也叫做初值问题(Initial Value Problem (IVP))。我们后面要讨论的数值方法都是围绕着 IVP 问题的。当然,常微分方程还有一类是边界问题(Boundary Value Problem(BVP)),这类问题我们留到偏微分方程的数值方法时一起讨论。

   那么,对于上面这样的一个微分方程,我们可以从两个不同的角度来得到几乎相同的数值方法。一种是将左边的导数进行泰勒展开的近似,另一种则是将整个方程积分。

1. 数值方法——泰勒展开观点

   我们考虑一个很小的值 $h$

  1. 对 $y(t+h)$ 做泰勒展开可以得到 $y(t+h)=y(t)+hy'(t)+\frac{h^2}{2!}y''(t)+\cdots$。我们舍去高次项,仅保留一阶导数,则 $y'(t)\approx \frac{y(t+h)-y(t)}{h}= f(y(t),t)$。这样就得到了Forward Euler方法,它的离散化误差(Discretization error)是 $e=\mathcal{O}(h)$。很显然,这个更新过程这是一个显式方法(Explicit method)。
  2. 同样的,对 $y(t)$ 做泰勒展开还可以得到 $y(t)=y(t+h)-hy'(t+h)+\frac{h^2}{2!}y''(t+h)-\cdots$。我们同样只保留一阶导数,就得到了 $y'(t+h)\approx \frac{y(t+h)-y(t)}{h}= f(y(t+h),t+h)$。这就是Backward Euler方法,它的离散化误差(Discretization error)同样是 $e=\mathcal{O}(h)$。但是,由于 $y(t+h)$ 同时出现在等式左右两端,因此需要求解这个方程,才能得到 $y(t+h)$ 的值,因此这是一个隐式方法(Implicit method)。
  3. 我们用 1)和 2)方法的平均值,就得到了一个新的方法:$ \frac{y(t+h)-y(t)}{h}= \frac{1}{2}\left(f(y(t+h),t+h)+f(y(t),t)\right)$。这个方法事实上就是所谓的梯形公式(Trapezoidal rule)。它仍然是一个隐式方法,但他的离散化误差提高到了是 $e=\mathcal{O}(h^2)$。有兴趣的小伙伴可以自己用泰勒展开验证一下。

2. 数值方法——数值积分观点

   上面所述的三种方法也可以从数值积分中推导出,我们先将微分方程两边同时积分,得到 $y(t+h)=y(t)+\int_t^{t+h}f \rm{d}s$。这个积分的值简单理解就是如下图所示,求函数 $f$ 与横坐标轴在 $t$ 到 $t+h$ 之间围成的面积。它们的数值方法有很多形式,例如

  1. 我们可以用下面蓝色的长方形区域来近似这个面积
    图
    图 1:用长方形近似面积
    $y(t+h)\approx y(t)+hf(y(t),t)$。这正好就是 Forward Euler 方法。
  2. 我们也可以用下面的绿色区域来近似这个面积
    图
    图 2:用长方形近似面积
    $y(t+h)\approx y(t)+hf(y(t+h),t+h)$。这就是Backward Euler方法。
  3. 还可以用下面的橙色梯形区域来近似
    图
    图 3:用梯形近似面积
    $y(t+h)\approx y(t)+\frac{1}{2}h\left(f(y(t),t)+f(y(t+h),t+h)\right)$。这正好就Trapezoidal rule,即梯形公式名字的由来。

3. 总结一下

   我们把 $y(t_i)$ 称作精确解或者解析解,它的数值近似解记做 $y_i$。

4. 求解步骤

   我们来简单概括一下常微分方程的求解过程:

  1. 从初始值开始,令 $y_0=y(t_0)=\hat{y}$。
  2. 根据所选用的数值方法,计算出 $y_1$。这里的方法可以是 Forward/Backward Euler,也可以是梯形公式,还可以是其他的一些方法。
  3. 将第 2 步的操作重复应用到由 $y_i$ 计算 $y_{i+1}$ 的过程。直到 $t_i\ge T$ 时,停止。

需要注意

   在上面的第 2 步使用 Backward Euler 或者梯形公式时,由于它们均为隐式方法,每一步都需要求解一个线性或非线性方程。我们考虑简单的常微分方程 $y'=f(y)$,以 Backward Euler 为例,数值方法为 $y_{n+1}=y_{n}+hf(y_{n+1})$ 即 $y_{n+1}-hf(y_{n+1})=y_{n}$

   这也就意味着,通常情况下如果采用隐式方法会比显式方法的运算复杂度更高。因此,在没有特殊要求的前提下,我们更偏向于使用显式方法求解常微分方程。(当然,我们会在后面的章节中具体讨论,哪些特殊情况下必须或者更偏向于使用隐式方法)

5. 更多方法

   也正是因为上面的原因,Karl Heun 将梯形公式改进成了显式方法----Heun 方法。它的核心思想就是用 Forward Euler 方法求出 $\tilde{y}_{n+1}=y_n+hf(y_n,t_n)$,然后带入到梯形公式的右边,$y_{n+1}=y_n+\frac{h}{2}\left(f(y_n,t_n)+f(\tilde{y}_{n+1},t_{n+1}) \right)$。这样整个过程的每一步都是显式方法,从而避免了任何可能出现的解方程的过程。当然,有些资料上面也会这样描述 Heun 方法:

\begin{equation} \begin{aligned} &k_1=f(y_n,t_n)~,\\ &k_2=f(y_n+hf(y_n,t_n),t_n+h))=f(y_n+hk_1,t_n+h)~,\\ &y_{n+1}=y_n+\frac{h}{2}(k_1+k_2)~. \end{aligned} \end{equation}

   事实上,这种形式是另一种更有名的方法的二阶特例,即龙格-库塔法(Runge-Kutta)。例如,我们最常见的 RK4 的形式为

\begin{equation} \begin{aligned} &k_1=f(y_n,t_n)~,\\ &k_2=f\left(y_n+\frac{h}{2}k_1, t_n+\frac{h}{2}\right)~,\\ &k_3=f\left(y_n+\frac{h}{2}k_2, t_n+\frac{h}{2}\right)~,\\ &k_4=f\left(y_n+hk_3, t_n+h\right)~,\\ &y_{n+1}=y_{n}+\frac{h}{6}(k_1+2k_2+2k_3+k_4)~. \end{aligned} \end{equation}

   当然,关于不同阶数 RK 方法的系数推导以及误差分析,已经超出了这个专题的范围,有兴趣的小伙伴可以留言,我会单独开一个专题来讨论。

最后一个例子

   我们来尝试着用上面的一些显式方法来求解一个非线性常微分方程:

\begin{equation} \begin{cases} y'(t) = y-\frac{1}{2}e^{\frac{t}{2}}\cdot \sin\left(5t\right) +5e^{\frac{t}{2}}\cdot \cos\left(5t\right) , 0\le t\le \pi\\ y(0)=0 \end{cases}~. \end{equation}

   我们分别用了 Forward Euler,Heun 和 RK4 方法,其中,Python 的科学计算包 scipy.integrate.solve_ivp 提供了例如 RK45(即我们上面提到的 RK4),BDF(后面我们会专门讨论它)以及其他的一些常用方法。

   有兴趣的小伙伴可以尝试着自己调整步长参数和方法,来看看误差的变化。

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np

def rhsODEs(t, y):
    return y - 0.5*np.exp(0.5*t)*np.sin(5*t)+5*np.exp(0.5*t)*np.cos(5*t)

# initial condition
y0 = [0]
N = 20

# Time steps
t_span = (0, np.pi)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# Solve for the ODE with R-K method
sol_ex = solve_ivp(rhsODEs, t_span, y0, method='RK45', t_eval=t_eval)
sol_fe = EulerForward(rhsODEs, t_span, y0, N)
sol_he = Heun(rhsODEs, t_span, y0, N)
t_evalRK = np.linspace(t_span[0], t_span[1], N)
sol_rk = solve_ivp(rhsODEs, t_span, y0, method='RK45', t_eval=t_evalRK)

# plot
fig, ax = plt.subplots(1, figsize=(6, 6))
ax.plot(sol_ex.t,sol_ex.y.T )
ax.plot(sol_fe[0], sol_fe[1],'-*' )
ax.plot(sol_he[0], sol_he[1],'-o' )
ax.plot(sol_rk.t,sol_rk.y.T, '-d')

ax.autoscale(enable=True, axis='both', tight=True)
ax.set_ylabel(r'$y(t)$')
ax.set_xlabel(r'$t$')
ax.legend(['Exact solution(RK45)','Forward Euler Method',
    'Heuns Method', 'Classical Runge-Kutta'])

   三种方法的结果如图,另外为了作为参照,我们选用了非常小的时间间隔,用 RK45 模拟出了一个解析解。

图
图 4:运行结果

   注:这里我自己写了 Forward Euler(或者叫 Euler Forward)和 Heun 方法如下:

import numpy as np

def EulerForward(func, t_span, y0, n):
    """ Explicit Euler method
    Parameters
    ----------
    func : function, function handler of the right hand side of the ODE(s);
    t_span : list, t_span[0] is the initial time, t_span[1] is the final time;
    y0 : float, list, initial condition(s) of the ODE(s);
    n : number of time steps, dt = (t_span[1]-t_span[0])/n.
        
    Returns
    -------
    t : list, the resulting time steps;
    y : list, the solutions of the ODE(s)
    """
    t = np.linspace(t_span[0], t_span[1], n)
    dt = t[1]-t[0]
    y = [y0]
    for i in range(n-1):
        ynp1 = y[i] + dt*func(t[i], y[i])
        y.append(ynp1)
    return t, y
    
    
def Heun(func, t_span, y0, n):
    """ Heun's method
    Parameters
    ----------
    func : function, function handler of the right hand side of the ODE(s);
    t_span : list, t_span[0] is the initial time, t_span[1] is the final time;
    y0 : float, list, initial condition(s) of the ODE(s);
    n : number of time steps, dt = (t_span[1]-t_span[0])/n.
        
    Returns
    -------
    t : list, the resulting time steps;
    y : list, the solutions of the ODE(s)
    """
    t = np.linspace(t_span[0], t_span[1], n)
    dt = t[1]-t[0]
    y = [y0]
    for i in range(n-1):
        k1 = func(t[i], y[i])
        k2 = func(t[i+1], y[i]+dt*k1)
        ynp1 = y[i] + 0.5*dt*(k1+k2)
        y.append(ynp1)
    return t, y

最后一点

   我们在这个例子中可以观察到一点,这三种方法的误差排序大体上是 Forward Euler > Heun > RK4。那么,下一章我们会来具体分析这其中的原因,也就是所谓的局部残差(local truncation error)和全局误差(global error)。


1. ^ 本文经作者同意转载自知乎专栏 《科学计算》,格式有少量修改。

                     

© 小时科技 保留一切权利