Lorenz 混沌系统的数值模拟

Lorenz 系统

很多微分方程都是确定性的：只要给定初始条件，理论上后面的状态也就确定了。 但“确定”并不等于“长期可预测”，还有著名的混沌系统，Lorenz 系统就是一个很经典的例子。

Lorenz 系统是三个变量组成的 ODE 系统：

$$ \begin{aligned} \dot x &= \sigma(y-x),\\ \dot y &= x(\rho-z)-y,\\ \dot z &= xy-\beta z. \end{aligned} $$

最常见的一组参数为 $\sigma=10, \rho=28, \beta=\frac{8}{3}$，初始条件取 $(x_0,y_0,z_0)=(1,1,1)$。

问题：一个完全确定的三维微分方程，为什么长期轨迹却很难预测？如果只是把初始值改动一点点，或者把数值积分的步长改小一点，会发生什么？

下面做几个简单的数值实验来观察这些现象。

轨迹可视化

首先直接对 Lorenz 方程使用 DOP853 进行求解。

DOP853：一种显式 RK8，搭配了 5 阶和 3 阶的误差估计以及自适应步长机制，适合求解刚性问题

得到的三维轨迹就是经典的 Lorenz attractor

其中的分量 $x(t)$ 如下

可以看到它的 $x$ 坐标一直在振荡，但是并没有一个很明显的固定周期。

由此可见，三维轨迹会在左右两个区域之间来回切换，看起来有点像蝴蝶的两片翅膀（“蝴蝶效应” 可能就是据此命名的？）。 虽然这里的运动轨迹很复杂，但轨迹却不会发散，长期被限制在一个有限区域内。

注意：

• 一条轨迹看起来复杂，只能说明运动比较复杂，并不能充分说明系统是混沌的。
• 在画吸引子时通常会丢掉前面一段时间：因为刚开始积分时，轨迹仍然明显受到人为给定初始条件的影响。过一段时间之后，轨迹才进入长期吸引子附近。 因此，本文画吸引子时去掉了前 10 个时间单位，这一段通常叫做 transient。

初值扰动的敏感性

我们考虑对初值进行小扰动，两组初值分别为

$$ (1,1,1), \quad (1+10^{-8},1,1). $$

初始状态的差别非常小，相对误差只有一亿分之一。

两组初值得到的轨迹的 $x(t)$ 如下图所示

观察可知：在前面一段时间里，两条曲线几乎完全重合，但是随着时间增加，它们会逐渐失去同步，最后变成完全不同的轨迹。

为了更清楚地对比，可以计算两个三维状态之间的距离：

$$ d(t)=\|u_1(t)-u_2(t)\|_2. $$

距离随时间的演化曲线如下图（对纵轴使用的是对数坐标）

由此可见，初值的微小扰动在一段时间后会快速扩大，经过一段时间后可能变成截然不同的两条轨迹，这就是初值敏感性。

对混沌系统的定性判断，其实还需要计算系统的三个 Lyapunov 指数：

• Lyapunov 指数大于 0，表示某个方向上的微小扰动平均会指数增长；
• 小于 0，表示扰动平均会衰减；
• 接近 0，表示既不明显增长也不明显衰减。

这里我们不对 Lyapunov 指数的定义和计算进行讨论，只是给出简单的说明：在当前系统的三个 Lyapunov 指数中，第一个指数大于 0，第二个非常接近 0，第三个小于 0。

数值求解的不可靠性

混沌系统的一个最重要的影响是：数值方法无法可靠地预测长期轨迹，因为数值方法的各种计算误差也相当于小扰动，可能会被混沌动力学不断放大。

在微分方程的数值求解中，这种情况叫做问题本身是不适定的，不适定性会给数值求解带来本质上的困难，通常我们只考虑适定问题的数值求解。

前面我们把初始条件加了小扰动，两条轨迹在长时间演化时会产生显著差异。现在我们继续说明，在减少数值求解的步长时，只能保证对应解轨迹的短期收敛性，长期演化后仍然会分离。

对于混沌问题的数值求解，虽然永远无法期望得到长期可靠的结果，但是至少高精度算法可以把明显失真的时间推得更远，而且在短时间内还是可以获得对应的数值精度的。

下面使用不同步长的精度 RK4 方法进行求解，比较不同步长得到的轨迹误差（与 DPO853 方法获取的参考解对比），可以看到，即使加密步长，在长期演化后仍然存在显著误差。

完整源码

import numpy as np
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

plt.rcParams.update(
    {
        "savefig.dpi": 300,
        "figure.autolayout": True,
        "text.usetex": True,
        "font.family": "serif",
        "font.serif": ["Times New Roman"],
        "pdf.fonttype": 42,
        "ps.fonttype": 42,
        "font.size": 14,
        "axes.labelsize": 14,
        "axes.titlesize": 14,
        "xtick.labelsize": 14,
        "ytick.labelsize": 14,
        "legend.fontsize": 14,
        "lines.linewidth": 1.0,
        "lines.markersize": 5,
        "axes.linewidth": 1.0,
        "lines.markerfacecolor": "none",
        "lines.markeredgecolor": "auto",
        "xtick.direction": "in",
        "ytick.direction": "in",
        "legend.facecolor": "white",
        "legend.edgecolor": "black",
        "legend.framealpha": 1.0,
    }
)

# ============================================================
# Configuration
# ============================================================

SIGMA = 10.0
RHO = 28.0
BETA = 8.0 / 3.0

U0 = np.array([1.0, 1.0, 1.0], dtype=float)

REF_METHOD = "DOP853"
REF_RTOL = 1e-11
REF_ATOL = 1e-13


def lorenz(t, u):
    x, y, z = u

    dx = SIGMA * (y - x)
    dy = x * (RHO - z) - y
    dz = x * y - BETA * z

    return np.array([dx, dy, dz], dtype=float)


def solve_reference(u0=U0, T=60.0, dt=0.01):
    t_eval = np.linspace(0.0, T, int(round(T / dt)) + 1)

    sol = solve_ivp(
        lorenz,
        (0.0, T),
        np.asarray(u0, dtype=float),
        t_eval=t_eval,
        method=REF_METHOD,
        rtol=REF_RTOL,
        atol=REF_ATOL,
    )

    if not sol.success:
        raise RuntimeError(sol.message)

    return sol.t, sol.y


# ============================================================
# Basic trajectory and attractor
# ============================================================


def basic_plots():
    t, u = solve_reference(T=60.0, dt=0.01)
    x, y, z = u

    # x(t)
    plt.figure(figsize=(10, 4))
    plt.plot(t, x, linewidth=0.8)
    plt.xlabel("$t$")
    plt.ylabel("$x(t)$")
    plt.title("Lorenz time series")
    plt.grid(alpha=0.25)
    plt.savefig("timeseries_x.png")
    plt.close()

    # Remove transient when plotting attractor
    mask = t >= 10.0

    fig = plt.figure(figsize=(8, 6))
    ax = fig.add_subplot(111, projection="3d")
    ax.plot(x[mask], y[mask], z[mask], linewidth=0.45)
    ax.set_xlabel("$x$")
    ax.set_ylabel("$y$")
    ax.set_zlabel("$z$")
    ax.set_title("Lorenz attractor")
    plt.savefig("lorenz_attractor.png")
    plt.close()


# ============================================================
# Sensitivity to initial conditions
# ============================================================


def sensitivity_experiment():
    perturbation = 1e-8

    u1 = U0.copy()
    u2 = U0.copy()
    u2[0] += perturbation

    U_init = np.concatenate([u1, u2])

    T = 50.0
    dt = 0.005
    t_eval = np.linspace(0.0, T, int(round(T / dt)) + 1)

    def twin_rhs(t, U):
        return np.concatenate(
            [
                lorenz(t, U[:3]),
                lorenz(t, U[3:]),
            ]
        )

    sol = solve_ivp(
        twin_rhs,
        (0.0, T),
        U_init,
        t_eval=t_eval,
        method=REF_METHOD,
        rtol=1e-11,
        atol=1e-13,
    )

    if not sol.success:
        raise RuntimeError(sol.message)

    u1_hist = sol.y[:3]
    u2_hist = sol.y[3:]
    distance = np.linalg.norm(u1_hist - u2_hist, axis=0)

    # Two x(t) curves
    plt.figure(figsize=(10, 4))
    plt.plot(sol.t, u1_hist[0], linewidth=0.8, label="$x_0 = 1$")
    plt.plot(
        sol.t,
        u2_hist[0],
        linewidth=0.8,
        label="$x_0 = 1 + 10^{-8}$",
    )
    plt.xlabel("$t$")
    plt.ylabel("$x(t)$")
    plt.title("Sensitivity to initial conditions")
    plt.legend()
    plt.savefig("initial_condition_x.png")
    plt.close()

    # Distance between the full 3D states
    plt.figure(figsize=(10, 4))
    plt.semilogy(
        sol.t,
        np.maximum(distance, 1e-16),
        linewidth=0.9,
    )
    plt.xlabel("$t$")
    plt.ylabel("$||u_1(t)-u_2(t)||_2$")
    plt.title("Growth of a $10^{-8}$ initial perturbation")
    plt.grid(alpha=0.25)
    plt.savefig("initial_condition_distance.png")
    plt.close()


# ============================================================
#  Long-time numerical divergence
# ============================================================


def rk4_integrate(f, u0, T, dt):
    def rk4_state_step(u, dt):
        k1 = f(0.0, u)
        k2 = f(0.0, u + 0.5 * dt * k1)
        k3 = f(0.0, u + 0.5 * dt * k2)
        k4 = f(0.0, u + dt * k3)

        return u + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)

    n = int(round(T / dt))
    dt = T / n

    t = np.linspace(0.0, T, n + 1)
    u = np.empty((3, n + 1), dtype=float)
    u[:, 0] = np.asarray(u0, dtype=float)

    for i in range(n):
        u[:, i + 1] = rk4_state_step(u[:, i], dt)

    return t, u


def numerical_divergence_experiment():
    T = 50.0
    dts = [0.02, 0.01, 0.005]

    ref = solve_ivp(
        lorenz,
        (0.0, T),
        U0,
        method=REF_METHOD,
        rtol=1e-13,
        atol=1e-15,
        dense_output=True,
        max_step=0.005,
    )

    horizons = []

    plt.figure(figsize=(10, 5))

    for dt in dts:
        t, u = rk4_integrate(lorenz, U0, T, dt)
        u_ref = ref.sol(t)

        error = np.linalg.norm(
            u - u_ref,
            axis=0,
        )

        plt.semilogy(
            t,
            np.maximum(error, 1e-16),
            linewidth=0.9,
            label=f"RK4 $h={dt:g}$",
        )

        idx = np.where(error > 1.0)[0]

        if idx.size:
            horizon = float(t[idx[0]])
        else:
            horizon = np.nan

        horizons.append((dt, horizon))

    plt.xlabel("$t$")
    plt.ylabel("distance from DOP853 reference")
    plt.title("Long-time numerical trajectory divergence")
    plt.legend()
    plt.grid(alpha=0.25)
    plt.savefig("long_time_numerical_divergence.png")
    plt.close()

    return horizons


basic_plots()

sensitivity_experiment()

numerical_divergence_experiment()