注解

此笔记本可在此处下载： 00_ODE.ipynb

代码作者：Ludovic Charleux <ludovic.charleux@univ-smb.fr>

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd
import numba
from scipy import integrate
from IPython.display import Image
%matplotlib notebook
plt.ioff()

前提条件和目标

必备条件

计算一阶或二阶线性可分离变量微分方程的精确解
集成
物理基本定律

目标

变换一个阶微分方程的步骤 \(p \in \mathbb{{N}}\) 转化为1阶系统
要识别颂歌的性质，表征它的状态变量
使用Euler，Taylor和Runge Kutta的方法
了解它们各自的优缺点，选择最适合所考虑问题的方法。
对它们进行编程并了解Python中最常用的预定义解算器。
分析和解释模型的结果

微分方程

广泛应用于物理
仅适用于特定情况的封闭式解决方案
需要数值求解器

它们与所研究的物理问题的数学公式相对应。要找到连续问题的精确解通常是很复杂的。数值方法的兴趣在于能够逼近问题的解。

我们将记住从物理问题到解决它的计算程序的两个关键步骤：

建模 (由ODE系统描述的数学模型)
数值逼近 (从可在计算机上编程的连续问题公式过渡到离散化问题公式)

两类微分方程

常微分方程

未知函数仅关于单个变量时间的导数 \(t\) 例如。

示例： 一维谐振子方程

\[\ dfrac d^2x dt^2+2\zeta\omega_0\dfrac dx dt+omega_0 x=0\]

这是一个常系数的二阶线性方程。

偏微分方程

未知函数关于几个变量、时间的导数 \(t\) 和空间 \((x, y, z)\) 例如。需要使用本课程中未介绍的专用方法(例如，有限元方法)

示例： 热方程

\[\ rho c_p\dfrac \部分t \部分t-k\delta t+s=0\]

这是一个抛物线方程(描述物体的热行为随时间的演变)。

微分方程是用来描述我们周围的世界的。它们在各个领域都可以找到，无论是机械、电子、经济等。。可以使用两种方法来解决这些问题：

这个 代数法 给出了精确解的表达式。然而，这种方法并不适用于所有类型的方程；
这个 数值方法 给出了一个柯西问题(ODE+I.C)解的近似解。

在下面，我们将集中讨论如何使用数值方法来求解常微分方程组。

引言示例：弹簧-质量系统

图像0

Source : Wikipedia

安装程序

点质量 \(m\) ，
弹簧刚度 \(k\) 。

问题表述：

\[\ ddot x+\omega_0^2 x=0\]

使用：

\[\欧米伽0=\sqrt \dfrac k m\]

封闭式解决方案

我们假设：

\[\left\lbrace\begin{align*}\]

那么解决方案是：

\[x（t）=x_0\cos（\omega_0 t）\]

def Solution(t):
    """
    Closed form solution.
    """
    return a0 * np.cos(omega0 * t)

tmax = 1.
a0 = 1.
omega0 = 2. * np.pi
ts = np.linspace(0., tmax, 200)
xs = Solution(ts)

plt.figure()
plt.plot(ts, xs, "k--", label = "Closed Form Solution")
plt.grid()
plt.xlabel("Time, $t$ [s]")
plt.ylabel("Amplitude, $x$ [m]")
plt.legend(loc = "best");

plt.show()

<IPython.core.display.Javascript object>

重新配方

让我们假设：

\[十=\]

因此：

\[\点X=\begin{bmatrix}\]

def Derivative(X, t):
    """
    ODE
    """
    return np.array([X[1], -omega0**2 * X[0]])

ODE的数值积分

利弊

优点：

任何颂歌都可以解决。
通用解决方案

缺点：

近似解
大多数情况下需要一台计算机

时间离散化

时间离散化 : \(t_0\) ， \(t_1\) ， \(\ldots\) ， \(t_i\) \(\ldots\) ， \(t_{{N-1}}\)
时间步长 \(\Delta t = t_{{i+1}} - t_i\) ，

Image(filename='time_grid.png')

显式欧拉法

直观的，
快地,
收敛缓慢，

\[X_{i+1}=X_i+\Δt\，f(X，t_i)\]

def Euler(func, X0, t):
    """
    Euler solver.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X  = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt-1):
        X[i+1] = X[i] + func(X[i], t[i]) * dt
    return X

Nt = 20
tv = np.linspace(0., tmax, Nt)
Xe = Euler(Derivative, [1., 0.], tv)
tl = ["$t_{{{0}}}$".format(i) for i in range(Nt)]

plt.figure()

plt.plot(ts, xs, "k--", label = "Closed Form Solution")
plt.plot(tv, Xe[:, 0], "or-", label = "Euler")
plt.grid()
plt.xlabel("Time, $t$ [s]")
plt.ylabel("Amplitude, $x$ [m]")
plt.legend(loc = "best")
plt.xticks(tv, tl)
plt.show()
#plt.close()

<IPython.core.display.Javascript object>

解的准确性

markers = "osv*<>"
Nt_values = np.array([5, 10, 20, 100])
plt.figure()
for i in range(len(Nt_values)):
    tvi = np.linspace(0., tmax, Nt_values[i])
    Xei = Euler(Derivative, [1., 0.], t =  tvi)
    xsi = Solution(tvi)
    #plt.plot(tvi, (Xei[:, 0] - xsi) / a0 * 100,
             #markers[i] + "-", label = "Euler, Nt = {0}".format(Nt_values[i]))
    plt.plot(tvi, np.abs(Xei[:, 0] - xsi), markers[i] + "-", label = "Euler, Nt = {0}".format(Nt_values[i]))
plt.grid()
plt.xlabel("Time, $t$ [s]")
plt.ylabel("Error")
plt.legend(loc = "best");

plt.show()
plt.close()

<IPython.core.display.Javascript object>

收敛顺序

\[E=C(\增量t)^p\]

因此，通过对数，我们得到以下等式：

\[\下括号{\log(E)}_{y}=C+p\下括号{\log(\Delta t)}_{x}\]

markers = "osv*<>"
Nt_values = np.logspace(1., 4., 20).astype(np.int32)
dt_values = tmax / (Nt_values -1)
Error = []
for i in range(len(Nt_values)):
    tvi = np.linspace(0., tmax, Nt_values[i])
    Xei = Euler(Derivative, [1., 0.], t =  tvi)
    xsi = Solution(tvi)
    e = (Xei[:, 0] - xsi).std()
    Error.append(e)

plt.figure()
plt.plot(dt_values, Error, "or-", label = "Explicit Euler error")
plt.plot(dt_values, dt_values, label = "$\Delta t$", linestyle ="--")
plt.grid()
plt.legend()
plt.xlabel("Times step, $\Delta t$ [s]")
plt.ylabel("Error (Standard Deviation)")
plt.yscale("log")
plt.xscale("log")

plt.show()
plt.close()

<IPython.core.display.Javascript object>

显式欧拉方法的阶是1

龙格库塔4

Wikipedia

具有以下特点的欧拉积分器的发展：

多点坡度评估（此处为4）
精心选择的权重与简单的解决方案相匹配。

\[X_{i+1}=X_i+\dfrac{\Delta t}{6}\Left(k_1+2k_2+2k_3+k_4\right)\]

使用：

\(k_1\) 是基于间隔开始处的斜率的增量，使用 \(X\) (欧拉法)；
\(k_2\) 是基于间隔中点处的斜率的增量，使用 \(\displaystyle X + \frac{{\Delta t}}{{2}} \times k_1\) ；
\(k_3\) 也是基于中点处的坡度的增量，但现在使用 \(\displaystyle X + \frac{{\Delta t}}{{2}} \times k_2\) ；
\(k_4\) 是基于间隔末尾的斜率的增量，使用 \(\displaystyle X + \Delta t \times k_3\) 。

def RK4(func, X0, t):
    """
    Runge Kutta 4 solver.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X  = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt-1):
        k1 = func(X[i], t[i])
        k2 = func(X[i] + dt/2. * k1, t[i] + dt/2.)
        k3 = func(X[i] + dt/2. * k2, t[i] + dt/2.)
        k4 = func(X[i] + dt    * k3, t[i] + dt)
        X[i+1] = X[i] + dt / 6. * (k1 + 2. * k2 + 2. * k3 + k4)
    return X

Xrk4 = RK4(Derivative, [1., 0.], tv)
plt.figure()
plt.plot(ts, xs, "k--", label = "Closed Form Solution", linewidth = 2.)
plt.plot(tv, Xe[:, 0], "or-", label = "Euler")
plt.plot(tv, Xrk4[:, 0], "sg-", label = "RK4")
plt.xticks(tv, tl)
plt.grid()
plt.xlabel("Time, $t$ [s]")
plt.ylabel("Amplitude, $x$ [m]")
plt.legend(loc = "best");

plt.show()
plt.close()

<IPython.core.display.Javascript object>

Odeint（带坐骨神经痛）

来自开源Fortran库的最新解决方案 ODEPACK .

它在科学类库有售 scipy.integrate.odeint .

参见 documentation

Xodeint = integrate.odeint(Derivative, [1., 0.], tv)
plt.figure()
plt.plot(ts, xs, "k--", label = "Closed Form Solution", linewidth = 2.)
plt.plot(tv, Xe[:, 0], "or-", label = "Euler")
plt.plot(tv, Xrk4[:, 0], "sg-", label = "RK4")
plt.plot(tv, Xodeint[:, 0], "*m-", label = "ODEint")
plt.xticks(tv, tl)
plt.grid()
plt.xlabel("Time, $t$ [s]")
plt.ylabel("Amplitude, $x$ [m]")
plt.legend(loc = "best");

plt.show()
plt.close()

<IPython.core.display.Javascript object>

请注意，在Scipy模块中还有其他用于求解ODE的原生Python方法，例如

scipy.integrate.solve_ivp or scipy.integrate.ode

隐式方法

动机

方程刚度
溶液的稳定性

后向欧拉

利弊

优点：

更准确

缺点：

必须解一个方程式
由于需要使用递归求根方法(割线法、牛顿-拉夫森法)，因此实现起来更加困难。

\[X_{i+1}=X_{i}+\Δt\，f(X，t_{i+1})\]

from scipy.optimize import fsolve

def EulerImp(func, X0, t):
    """
    Euler Implicit solver.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X  = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt-1):
        func_eq = lambda Xnp : Xnp - X[i] + dt * func(Xnp, t[i+1])
        X[i+1] = fsolve(func_eq, X[i])

    return X

Xeimp = EulerImp(Derivative, [1., 0.], tv)
plt.figure()
plt.plot(ts, xs, "k--", label = "Closed form solution", linewidth = 2.)
plt.plot(tv, Xe[:, 0], "or-", label = "Euler explicit")
plt.plot(tv, Xrk4[:, 0], "sg-", label = "RK4")
plt.plot(tv, Xodeint[:, 0], "*m-", label = "ODEint")
plt.plot(tv, Xeimp[:,0], ">b-", label="Euler implicit")
plt.xticks(tv, tl)
plt.grid()
plt.xlabel("Time, $t$ [s]")
plt.ylabel("Amplitude, $x$ [m]")
plt.legend(loc = "best");

plt.show()

<IPython.core.display.Javascript object>

Nt_values = np.logspace(1., 4., 20).astype(np.int32)
dt_values = tmax / (Nt_values -1)
ErrorE = []
ErrorRK4 = []
ErrorODEint = []
ErrorEImp = []
for i in range(len(Nt_values)):
    tvi = np.linspace(0., tmax, Nt_values[i])
    Xei = Euler(Derivative, [1., 0.], t =  tvi)
    Xrk4i = RK4(Derivative, [1., 0.], t =  tvi)
    Xodeinti = integrate.odeint(Derivative, [1., 0.], t =  tvi)
    XeImpi = EulerImp(Derivative, [1., 0.], t =  tvi)
    xsi = Solution(tvi)
    ErrorE.append((Xei[:, 0] - xsi).std())
    ErrorRK4.append((Xrk4i[:, 0] - xsi).std())
    ErrorODEint.append((Xodeinti[:, 0] - xsi).std())
    ErrorEImp.append((XeImpi[:, 0] - xsi).std())

plt.figure()
plt.plot(dt_values, ErrorE, "or-", label = "Euler")
plt.plot(dt_values, ErrorRK4, "sg-", label = "RK4")
plt.plot(dt_values, ErrorODEint, "*m-", label = "ODEint")
plt.plot(dt_values, ErrorEImp, "+b-", label = "EulerImp")
plt.plot(dt_values, dt_values, label = "$\Delta t$", linestyle ="--")
plt.plot(dt_values, dt_values**4, label = "$\Delta t^4$", linestyle ="--")
plt.grid()
plt.xlabel("Times step, $\Delta t$ [s]")
plt.ylabel("Error (Standard Deviation)")
plt.legend(loc = "best")
plt.yscale("log")
plt.xscale("log")

plt.show()

<IPython.core.display.Javascript object>

在本课程中，我们只看到了一步数值格式，即计算量的方法。 \(X_i\) 预计将是以下内容的近似值 \(X(t_i)\) 由以下人员定义：

\[\begin{split}X_{i+1}=X_{i}+\Δt\\φ(X，t_i，\Δt)\quad\for all\0\leq i\leq N-1\end{split}\]

\(\Phi(X, t_i, \Delta t) := f(X, t_i)\) 对于显式欧拉方法
\(\Phi(X, t_i, \Delta t) := f(X, t_{{i+1}})\) 隐式欧拉法
\(\Phi(X, t_i, \Delta t) := \displaystyle \frac{{1}}{{6}} \left(k_1 + 2 k_2 + 2 k_3 + k_4 \right)\) 对于RK4方法

对于那些想走得更远的人来说，还有另一类数值方法，那就是 multi-step methods 。启动一段时间后，可以使用以下近似值 \(X\) : \(X_{{i-1}}, X_{{i-2}},..., X_{{i-k+1}}\) 在前面的步骤中 \(X_i\) 。

下面是一些经典的例子：

亚当斯方法：亚当斯-巴什福斯、亚当斯-莫尔顿
Nyström
米尔恩-辛普森
向后微分公式(BDF)

Multistep methods