程序设计与计算思维

Computer Programming and Computational Thinking

第 18 讲：天气预报难题

2025—2026学年度春季学期

清华大学 韩文弢

问题引入：天气与气候

明天会下雨吗？我们应该能够利用对大气当前状态的观测，以及以偏微分方程形式表示的数学模型来预测这一点——这些模型可以预测风如何移动空气、温度将是多少、以及有多少水会凝结降水。

那么预测下周的天气呢？明年呢？

我们从经验中知道，我们只能在短期内预测天气，大约一周左右（而且仍然只有一定的概率）。甚至一个月后的预测似乎都是不可能的。

另一方面，我们可以讨论气候及其在更长时间尺度上的行为。例如，"北京四月份的平均气温是15.1°C"。短期天气具有显著的随机（或看似随机的）成分，但当我们观察平均值和标准差等统计性质时，就变得更具可预测性；这就是所谓的气候。

非线性动力学：稳定性与分岔

本讲将学习最简单的天气模型——由3个耦合的非线性常微分方程组成的洛伦兹方程组——可以表现出不可预测的行为，即混沌。

首先通过研究两个更简单的动力学模型（一维和二维）来逐步深入理解。

准备工作

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
import pandas as pd
from ipywidgets import interact, FloatSlider

plt.rcParams['font.family'] = 'Noto Sans CJK SC'

回顾：常微分方程（ODE）

回想一下，我们正在研究微分方程，目标是模拟气候随时间的演化。

最简单的此类模型是常微分方程（ODE），其中一个或少数几个连续变量随时间连续演化，模型指定了它们作为当前值函数的瞬时变化率。

不显式依赖于时间的ODE（即自治系统）的一般形式为

$$\frac{dx(t)}{dt} = f(x(t)),$$

或

$$\dot{x}(t) = f(x(t)),$$

带有初始条件 $x(t=0) = x_0$。

这里 $\dot{x}(t)$ 表示函数 $t \mapsto x(t)$ 在时刻 $t$ 的导数。

解析解与数值解

对于某些特殊形式的ODE，我们可以求出解析解（即用公式精确表达的解）：

方程	解析解
$\dot{x} = a$	$x(t) = x_0 + at$
$\dot{x} = a - bx$	$x(t) = \frac{a}{b} + \left(x_0 - \frac{a}{b}\right)e^{-bt}$
$\dot{x} = a - bx + f(t)$	积分因子法，一般无简洁形式

然而，大多数非线性ODE（如我们即将看到的 $\dot{x} = x(1-x)$）没有解析解，必须用数值方法近似求解。

t = np.linspace(0, 5, 200)
x0 = 1.0

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

a = 2.0
axes[0].plot(t, x0 + a * t, 'b-', lw=2)
axes[0].set_title(r'$\dot{x} = a$ (常数增长)')
axes[0].set_xlabel('$t$'); axes[0].set_ylabel('$x(t)$')

a, b = 3.0, 1.5
x_eq = a / b
axes[1].plot(t, x_eq + (x0 - x_eq) * np.exp(-b * t), 'r-', lw=2)
axes[1].axhline(y=x_eq, color='gray', ls='--', alpha=0.5, label=f'平衡态 $x^* = {x_eq:.1f}$')
axes[1].set_title(r'$\dot{x} = a - bx$ (指数趋近)')
axes[1].set_xlabel('$t$'); axes[1].legend()

a, b = 2.0, 1.0
sol = solve_ivp(lambda t, x: a - b * x + np.sin(t), [0, 5], [x0], t_eval=t)
axes[2].plot(sol.t, sol.y[0], 'g-', lw=2)
axes[2].set_title(r'$\dot{x} = a - bx + \sin(t)$ (含外力)')
axes[2].set_xlabel('$t$')

plt.tight_layout()
plt.show()

求解此类方程最简单的数值方法是欧拉法，在该方法中，我们将方程转化为显式的时间步进格式：取长度为 $h$ 的小时间步，将近似导数表示为

$$\frac{dx(t)}{dt} \simeq \frac{x(t + h) - x(t)}{h},$$

得到数值方法

$$x_{n+1} = x_{n} + h \, f(x_n),$$

其中 $x_n$ 是在第 $n$ 个时间步 $t_n$ 处对真实解 $x(t_n)$ 的近似。

一维：细菌增长模型

让我们用这个来模拟描述细菌种群动力学的简单非线性ODE。细菌以速率 $\lambda$ 繁殖，前提是有充足的食物，在这种情况下我们有 $\dot{x} = \lambda x$。但对充足食物的需求会将可持续种群限制在一个值 $K$，有时称为环境容纳量。

增长和饱和的综合效应的最简单模型如下：

$$\dot{x} = \lambda \, x \, (K - x).$$

当 $x$ 接近 $0$ 时，增长率为 $\lambda$，但随着 $x$ 的增加，增长率降低。

（这有时被称为逻辑斯谛微分方程。）

我们的目标是运用计算思维，但实际上我们并不太关心精确的时间动力学，而是关注系统行为的定性特征。例如，在长时间（形式上是 $t \to \infty$）后，种群是否会无限增大？或者它是否会（例如）在某个特定值附近振荡？或者收敛到某个特定大小？这就构成了非线性动力学或动力系统理论的研究主题。

使用欧拉法模拟系统来尝试猜测这个问题的答案。在实践中不应该使用欧拉法，而应该使用经过验证的库和提供更高精度的算法。

将变量重新缩放到最简单的形式是有用的：

$$\dot{x} = x \, (1 - x).$$

可以通过定义新的无量纲空间和时间变量 $x'$ 和 $t'$ 来实现：$x = K \, x'$，得到 $K \dot{x'} = \lambda K \, x' (K - K x')$，然后令 $t' = \lambda K \, t$，得到 $dx' / dt' = x' (1 - x')$。

定义一个表示ODE右端的函数：

def logistic(x):
    return x * (1 - x)

用欧拉法模拟并绘制轨迹 $x(t)$ 随时间 $t$ 的变化：

def euler_step(f, x, h):
    return x + h * f(x)

def euler(f, x0, h, t_final):
    ts = [0.0]
    xs = [x0]
    x = x0
    t = 0.0
    while t < t_final:
        x = euler_step(f, x, h)
        t += h
        xs.append(x)
        ts.append(t)
    return np.array(ts), np.array(xs)

ts, xs = euler(logistic, 0.5, 0.01, 20.0)

plt.figure(figsize=(6, 4))
plt.plot(ts, xs, lw=3)
plt.scatter([ts[0]], [xs[0]], zorder=5)
plt.xlabel('$t$')
plt.ylabel('$x(t)$')
plt.ylim(0.4, 1.1)
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(8, 5))

ts_q, xs_q = euler(logistic, 0.5, 0.01, 10.0)
ax.plot(ts_q, xs_q, 'b-', lw=3, label='$x(t)$')

T, X = np.meshgrid(np.linspace(0, 10, 20), np.linspace(0, 1.5, 12))
DT = np.ones_like(T)
DX = logistic(X)
M = np.sqrt(DT**2 + DX**2)

ax.quiver(T, X, DT/M, DX/M, M, cmap='coolwarm', alpha=0.5, scale=25)
ax.set_xlabel('$t$')
ax.set_ylabel('$x$')
ax.set_title(r'$\dot{x} = x(1-x)$ 的方向场与轨迹')
ax.set_xlim(0, 10)
ax.set_ylim(0, 1.5)
ax.legend()
plt.tight_layout()
plt.show()

上图中，每个箭头表示在该点 $(t, x)$ 处轨迹的斜率 $\dot{x}$。蓝色轨迹沿着箭头方向前进，最终趋向稳定不动点 $x^* = 1$。箭头颜色表示斜率大小：红色区域斜率大（变化快），蓝色区域斜率小（变化慢）。

我们看到对于这个特定的初始条件，解在一段时间后似乎稳定在一个固定值，然后保持在该值。

这样的值被称为ODE的不动点、驻点或稳态。

定性行为：不动点及其稳定性

x0 = 0.5

ts, xs = euler(logistic, x0, 0.01, 5.0)

plt.figure(figsize=(6, 4))
plt.plot(ts, xs, lw=3)
plt.scatter([0.0], [x0], zorder=5)
plt.axhline(y=0, linestyle='--', color='gray')
plt.xlabel('$t$')
plt.ylabel('$x(t)$')
plt.ylim(-1, 2)
plt.tight_layout()
plt.show()

让我们看看其他初始条件会发生什么：

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

for x0 in np.arange(-0.5, 2.1, 0.1):
    ts, xs = euler(logistic, x0, 0.05, 5.0)
    plt.plot(ts, xs, alpha=0.8, lw=1)

plt.xlabel('$t$')
plt.ylabel('$x(t)$')
plt.ylim(-1, 2)
plt.title('$\\dot{x} = x(1-x)$ 的轨迹')
plt.tight_layout()
plt.show()

/var/folders/1m/86zn99zs5p3fvt1c70k614n80000gn/T/ipykernel_88145/1817473811.py:2: RuntimeWarning: overflow encountered in scalar multiply
  return x * (1 - x)

我们看到所有从 $x_0=1.0$ 附近出发的曲线似乎都在长时间后收敛到1。如果系统精确地从0开始，它将永远停留在那里。然而，如果它从0附近开始（无论哪一侧），它都会远离0（在0的同一侧）——从负值开始，$x$ 会变得越来越负。（尽管负种群在原始的种群动力学解释中没有意义，但我们仍然可以研究具有负初始条件的方程动力学，因为它也可能模拟其他系统。）

特殊值 $x^*_1=1$ 和 $x^*_2=0$ 被称为微分方程的驻点或不动点。如果我们从 $x^*_i$ 开始，那里的导数为 $f'(x^*_i) = 0$，因此我们不能离开 $x^*_i$！不动点可以作为函数 $f$ 的零点或根找到，即满足 $f(x^*) = 0$ 的值 $x^*$。

我们看到，两种类型的不动点是定性不同的：从 $x^*_1 = 1$ 附近出发的轨迹会朝向 $x^*_1$ 移动，而从 $x^*_2 = 0$ 附近出发的轨迹会远离它移动。我们说 $x^*_1$ 是稳定不动点，$x^*_2$ 是不稳定不动点。

通常不可能找到不动点位置和稳定性的解析公式；相反，我们可以使用数值寻根算法，例如牛顿法。

状态空间：向量场与相图

Python 语法：try/except 异常处理

程序运行时可能遇到意外情况（如除以零、文件不存在等），Python 会抛出异常（exception）。如果不处理，程序会崩溃。

try/except 语句让我们捕获并处理异常，使程序能够优雅地继续运行：

try:
    result = risky_operation()   # 可能出错的代码
except SomeError as e:           # 捕获特定类型的异常
    print(f'出错了: {e}')         # 处理异常

可以同时捕获多种异常类型。如果不确定异常类型，可以用裸 except（但不推荐，因为会隐藏所有错误）。

Python 语法：三元表达式

在下面的代码中，需要根据条件选择不同的值。Python 提供了三元表达式（又称条件表达式），可以在一行中完成简单的条件判断：

值1 if 条件 else 值2

如果条件为真，表达式的结果为 值1，否则为 值2。这等价于：

if 条件:
    result = 值1
else:
    result = 值2

三元表达式让代码更简洁，常用于赋值、函数参数和列表推导式中。

如果我们想找到给定初始条件的完整轨迹，那么我们需要求解方程，无论是数值求解还是解析求解。

然而，我们可能只需要关于系统的较少信息，例如长时间或渐近动力学。事实证明，我们可以不显式求解ODE就获得一些关于它的信息！这就是研究非线性系统的定性方法。

与上面绘制轨迹 $x(t)$ 作为时间 $t$ 的函数不同，让我们使用一种不同的图形表示，绘制状态空间或相空间：这是所有可能因变量值（"状态"）的集合（"空间"）。对于上面的ODE，只有一个因变量 $x$，所以状态空间是实数轴 $\mathbb{R}$。

在 $x$ 的每个可能值处，ODE给出了该点处 $x(t)$ 变化率的信息。让我们在该点画一个箭头，指向放在该点的粒子会移动的方向：如果 $\dot{x} > 0$ 则向右，如果 $\dot{x} < 0$ 则向左。

def numerical_derivative(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

def find_zeros(f, x_min, x_max, n_points=10000):
    roots = []
    xs = np.linspace(x_min, x_max, n_points)
    vals = np.array([f(x) for x in xs])
    for i in range(len(vals) - 1):
        if vals[i] * vals[i+1] < 0:
            try:
                root = brentq(f, xs[i], xs[i+1])
                roots.append(root)
            except ValueError:
                pass  # brentq 未找到根时跳过
    return roots

def plot_vector_field_1d(f, x_range=(-2, 2.5), ax=None):
    if ax is None:
        fig, ax = plt.subplots(figsize=(8, 2))
    
    xlo, xhi = x_range
    arrow_size = 0.07
    tol = 1e-5
    
    for x in np.arange(xlo, xhi, 0.2):
        d = f(x)
        if d > tol:
            ax.annotate('', xy=(x + arrow_size, 0), xytext=(x - arrow_size, 0),
                        arrowprops=dict(arrowstyle='->', color='blue', alpha=0.5, lw=1.5))
        elif d < -tol:
            ax.annotate('', xy=(x - arrow_size, 0), xytext=(x + arrow_size, 0),
                        arrowprops=dict(arrowstyle='->', color='red', alpha=0.5, lw=1.5))
    
    roots = find_zeros(f, -10, 10)
    for root in roots:
        stability = numerical_derivative(f, root)
        color = 'green' if stability > 0 else 'black'
        marker = 's' if stability > 0 else 'o'
        size = 60 if stability > 0 else 80
        ax.scatter([root], [0], c=color, marker=marker, s=size, alpha=0.7, zorder=5)
    
    ax.set_xlim(xlo, xhi)
    ax.set_ylim(-0.5, 0.5)
    ax.set_yticks([])
    return ax

fig, ax = plt.subplots(figsize=(8, 2))
plot_vector_field_1d(logistic, ax=ax)
plt.title('logistic方程的向量场（绿色方块=不稳定，黑色圆点=稳定）')
plt.tight_layout()
plt.show()

这个向量场确实给出了动力学的定性图像。它没有告诉我们每个区域中动力学发生得多快，但它指示了趋势。我们根据稳定性对不动点进行了编码；这可以使用在不动点处计算的导数 $f'(x^*)$ 来计算，因为该导数控制了附近初始条件 $x^* + \delta x$ 的行为。

不稳定不动点显示为绿色方块，稳定不动点显示为灰色圆点；我们将在整个笔记本中使用这一约定。

分岔

现在假设系统中有一个可以变化的参数 $\mu$。对于 $\mu$ 的每个值，我们有一个不同的ODE

$$\dot{x} = f_\mu(x).$$

例如， $$\dot{x} = \mu + x^2.$$

让我们为 $\mu$ 的不同值绘制状态空间：

Python 语法：lambda 与闭包变量捕获

在下面的代码中，需要在循环中为每个 $\mu$ 值创建不同的函数传给 plot_vector_field_1d。这时 lambda（匿名函数）非常方便：

lambda 参数: 表达式

等价于用 def 定义一个只有 return 语句的函数，但更简洁。

但有一个常见的陷阱！ 在循环中创建 lambda 时，它捕获的是变量的引用而非值。循环结束后所有 lambda 都会使用同一个（最后的）值：

# 错误写法：所有 lambda 都输出 4
funcs = [lambda: i for i in range(5)]
print([f() for f in funcs])  # [4, 4, 4, 4, 4]

# 正确写法：用默认参数绑定当前值
funcs = [lambda i=i: i for i in range(5)]
print([f() for f in funcs])  # [0, 1, 2, 3, 4]

默认参数 i=i 在定义时求值，因此每个 lambda 都"记住"了当时的 i。在下面的代码中，lambda x, m=mu_val: g(m, x) 就用了这个技巧。

def g(mu, x):
    return mu + x**2

fig, axes = plt.subplots(1, 5, figsize=(16, 2))
mu_values = [-1.0, -0.5, 0.0, 0.5, 1.0]

for ax, mu_val in zip(axes, mu_values):
    plot_vector_field_1d(lambda x, m=mu_val: g(m, x), x_range=(-2, 2), ax=ax)
    ax.set_title(f'$\\mu = {mu_val}$')

plt.tight_layout()
plt.show()

现在让我们将所有向量场收集到一张图中。水平轴代表参数 $\mu$ 的不同可能值：

def plot_bifurcation_vector_field(f, mu_range=(-2, 2), x_range=(-2, 2)):
    fig, ax = plt.subplots(figsize=(8, 6))
    
    xlo, xhi = x_range
    arrow_size = 0.07
    tol = 1e-5
    
    for mu_val in np.arange(mu_range[0], mu_range[1], 0.15):
        for x in np.arange(xlo, xhi, 0.2):
            d = f(mu_val, x)
            if d > tol:
                ax.annotate('', xy=(mu_val, x + arrow_size), xytext=(mu_val, x - arrow_size),
                            arrowprops=dict(arrowstyle='->', color='blue', alpha=0.4, lw=1.2))
            elif d < -tol:
                ax.annotate('', xy=(mu_val, x - arrow_size), xytext=(mu_val, x + arrow_size),
                            arrowprops=dict(arrowstyle='->', color='red', alpha=0.4, lw=1.2))
        
        roots = find_zeros(lambda x: f(mu_val, x), -10, 10)
        for root in roots:
            stability = numerical_derivative(lambda x: f(mu_val, x), root)
            color = 'green' if stability > 0 else 'black'
            marker = 's' if stability > 0 else 'o'
            size = 15 if stability > 0 else 20
            ax.scatter([mu_val], [root], c=color, marker=marker, s=size, alpha=0.5)
    
    ax.set_xlim(*mu_range)
    ax.set_xlabel('$\\mu$')
    ax.set_ylabel('给定 $\\mu$ 下的不动点与动力学')
    ax.set_aspect('equal')
    return ax

plot_bifurcation_vector_field(g)
plt.title('$\\dot{x} = \\mu + x^2$ 的分岔图')
plt.tight_layout()
plt.show()

我们看到在临界值 $\mu_c = 0$ 处，系统的行为发生了定性变化：当 $\mu_c < 0$ 时有两个不动点，而当 $\mu_c > 0$ 时完全没有不动点。这种定性变化被称为分岔。在这种特定情况下，两个不动点在鞍-结或折叠分岔中碰撞。

一维系统：双稳态与滞后

现在让我们看看以下系统的动力学：

$$\dot{x} = \mu + x - x^3.$$

def h(mu, x):
    return mu + x - x**3

让我们绘制分岔图：

plot_bifurcation_vector_field(h)
plt.title('$\\dot{x} = \\mu + x - x^3$ 的分岔图')
plt.tight_layout()
plt.show()

我们看到在 $\mu$ 的某个取值范围内，存在三个共存的不动点，其中两个稳定，一个不稳定。由于系统可以停留在两个稳定不动点中，我们说系统是双稳态的。

现在我们理解了图的含义和动力学，让我们只绘制不动点 $x^*(\mu)$ 作为 $\mu$ 的函数。这样的图被称为分岔图：

def plot_fixed_points(f, mu_range=(-2, 2), x_range=(-2, 2)):
    fig, ax = plt.subplots(figsize=(8, 6))
    
    mu_values = np.arange(mu_range[0], mu_range[1], 0.01)
    
    for mu_val in mu_values:
        roots = find_zeros(lambda x: f(mu_val, x), x_range[0], x_range[1])
        for root in roots:
            stability = numerical_derivative(lambda x: f(mu_val, x), root)
            color = 'green' if stability > 0 else 'black'
            marker = 's' if stability > 0 else 'o'
            size = 10 if stability > 0 else 15
            ax.scatter([mu_val], [root], c=color, marker=marker, s=size, alpha=0.5)
    
    ax.set_xlabel('$\\mu$')
    ax.set_ylabel('$x^*(\\mu)$')
    ax.set_aspect('equal')
    return ax

plot_fixed_points(h)
plt.title('$\\dot{x} = \\mu + x - x^3$ 的不动点')
plt.tight_layout()
plt.show()

曲线的各段被称为分支。

滞后效应

假设我们现在考虑缓慢改变参数 $\mu$。如果我们稍微改变参数 $\mu$，系统不再处于不动点，因为当 $\mu$ 改变时不动点的位置会移动。然而，系统随后会弛豫：它将遵循新 $\mu$ 值下的动力学，并迅速收敛到该新 $\mu$ 值附近的新不动点。

例如，从 $\mu=-2$ 开始，系统将保持在下面的黑色（稳定）分支上，直到 $\mu \approx 0.4$。在那一点，两个不动点碰撞并湮灭！之后附近不再有不动点。然而，在更远处有另一个不动点，现在将吸引所有轨迹，因此系统快速过渡到那个不动点。

现在假设我们再次降低参数。系统现在将沿着上分支返回，直到 $\mu \approx -0.4$，此时它将再次跳回下面。

对于参数值 $\mu$ 在区间 $[-0.4, 0.4]$ 中的每个值，都存在双稳态，即两个不动点与同一 $\mu$ 值共存（连同第三个不可观测的不稳定不动点）。

系统追踪不同的稳定分支取决于我们从哪里开始，即取决于动力学的历史，这一事实被称为滞后效应。

这种滞后行为在许多科学和工程背景中都可以找到，包括生物学中的开关，例如基因开关，以及地球气候的历史动力学。

慢-快系统

当我们让参数 $\mu$ 变化时，我们实际上在做什么？实际上，我们现在有一个包含两个耦合ODE的系统，例如

$$\dot{x} = \mu + x - x^3;$$ $$\dot{\mu} = \epsilon,$$

其中 $\mu$ 以慢速 $\epsilon$ 变化。在远短于 $1/\epsilon$ 的时间尺度上，$x$ 的动力学"不知道"$\mu$ 在变化，因此它会收敛到当前 $\mu$ 值的不动点 $x^*(\mu)$。

然而，实际上 $\mu$ 确实逐渐变化，所以 $x$ 的值将有效地沿着曲线 $x(t) \simeq x^*(\mu(t))$ "滑动"，追踪 $\mu$ 变化时的不动点曲线。

一旦 $\mu$ 达到临界值 $\mu_c$，附近就不再有不动点可以移动，向量场将使系统快速过渡到远处的替代不动点。

如果我们现在反转 $\mu$ 的动力学，我们将沿着上分支滑回。

二维：化学反应中的振荡——布鲁塞尔子模型

分岔不局限于一维系统；事实上，高维系统中的动力学可以更加丰富。

例如，让我们看看布鲁塞尔子模型。这个模型模拟了一种振荡的化学反应，被称为化学钟：

$$\begin{aligned} \dot{x} &= a + x^2 y - bx - x \\ \dot{y} &= b x - x^2 y \end{aligned}$$

def brusselator(t, state, a, b):
    x, y = state
    dx = a + x**2 * y - b * x - x
    dy = b * x - x**2 * y
    return [dx, dy]

a_val = 1.0
b_val = 1.5

u0 = [1, 1]
t_span = (0, 50)
sol = solve_ivp(brusselator, t_span, u0, args=(a_val, b_val),
                max_step=0.01, dense_output=True)

print(f"Brusselator求解完成: t从{sol.t[0]:.1f}到{sol.t[-1]:.1f}, 共{len(sol.t)}个时间步")

Brusselator求解完成: t从0.0到50.0, 共5002个时间步

a_val = 1.0
b_val = 1.5
t_span = (0, 10)

fig, ax = plt.subplots(figsize=(7, 7))
ax.set_facecolor('black')

for x0 in np.arange(0, 5.5, 1.0):
    for y0 in np.arange(0, 5.5, 1.0):
        sol = solve_ivp(brusselator, t_span, [x0, y0], args=(a_val, b_val),
                        max_step=0.02, dense_output=True)
        ax.plot(sol.y[0], sol.y[1], lw=1.5, alpha=0.8)

# 绘制方向场
xs_arr = []
ys_arr = []
for x in np.arange(0, 5, 0.25):
    for y in np.arange(0, 5, 0.25):
        v = np.array(brusselator(0, [x, y], a_val, b_val))
        v = v / (np.linalg.norm(v) * 30)
        xs_arr.extend([x - v[0], x + v[0], np.nan])
        ys_arr.extend([y - v[1], y + v[1], np.nan])

ax.plot(xs_arr, ys_arr, alpha=0.5, color='gray', lw=0.8)

ax.set_xlim(0, 5)
ax.set_ylim(0, 5)
ax.set_aspect('equal')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_title('布鲁塞尔子模型的状态空间轨迹')
plt.tight_layout()
plt.show()

/var/folders/1m/86zn99zs5p3fvt1c70k614n80000gn/T/ipykernel_88145/2115798127.py:20: RuntimeWarning: invalid value encountered in divide
  v = v / (np.linalg.norm(v) * 30)

我们看到在 $b$ 的某个临界值以上，不动点变得不稳定并产生一个吸引的周期轨道。这就是Hopf分岔。

三维：洛伦兹方程中的混沌

洛伦兹方程构成了代表大气的流体层中对流的（非常）简化模型，由Edward Lorenz在1963年发表的一篇著名论文中首次研究，工作在MIT完成。它们代表了对混沌行为的开创性数值研究。

洛伦兹方程是一组三个耦合的ODE。它们具有看似简单的形式，只有两个非线性项：

$$\begin{align} \dot{x} &= \sigma (y - x) \\[6pt] \dot{y} &= x (\rho - z) - y \\[6pt] \dot{z} &= x y - \beta z \end{align}$$

尽管如此，我们将看到它们可以表现出极其复杂的行为。实际上，大多数具有至少三个变量的非线性ODE都可以预期表现出类似的复杂行为。

让我们使用 scipy.integrate.solve_ivp 求解该系统。我们将取经典参数值 $\sigma = 10$ 和 $\beta = 8/3$，但允许 $\rho$ 变化；其经典值为 $28$。

def lorenz(t, state, sigma, rho, beta):
    x, y, z = state
    dx = sigma * (y - x)
    dy = x * (rho - z) - y
    dz = x * y - beta * z
    return [dx, dy, dz]

sigma = 10.0
rho = 10.0
beta = 8.0 / 3.0

lorenz_sol = solve_ivp(lorenz, (0, 100), [0.01, 0.01, 0.01],
                       args=(sigma, rho, beta),
                       max_step=0.01, dense_output=True)

print(f"洛伦兹方程求解完成: rho={rho}, 共{len(lorenz_sol.t)}个时间步")

洛伦兹方程求解完成: rho=10.0, 共10001个时间步

rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0

lorenz_sol = solve_ivp(lorenz, (0, 100), [0.01, 0.01, 0.01],
                       args=(sigma, rho, beta),
                       max_step=0.01, dense_output=True)

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot(lorenz_sol.y[0], lorenz_sol.y[1], lorenz_sol.y[2],
        lw=0.5, alpha=0.8)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
ax.set_xlim(-25, 25)
ax.set_ylim(-25, 25)
ax.set_zlim(0, 60)
ax.set_title(f'洛伦兹吸引子 ($\\rho = {rho}$)')
plt.tight_layout()
plt.show()

随着 $\rho$ 的增加，我们看到一系列分岔。在临界值以上，轨迹收敛到一个分形的奇异吸引子，其上的动力学是混沌的。

确定性混沌发生在附近初始条件在状态空间中以指数速度分离的情况。这被称为蝴蝶效应：蝴蝶拍动翅膀对大气状态的扰动最终可能被放大到改变龙卷风移动的方向。

我们可以通过稍微扰动初始条件并计算两个系统之间的距离作为 $t$ 的函数来简单地看到这一点：

rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0
eps = 1e-10

lorenz_sol1 = solve_ivp(lorenz, (0, 50), [0.01, 0.01, 0.01],
                        args=(sigma, rho, beta),
                        max_step=0.01, dense_output=True)

lorenz_sol2 = solve_ivp(lorenz, (0, 50), [0.01 + eps, 0.01 + eps, 0.01 + eps],
                        args=(sigma, rho, beta),
                        max_step=0.01, dense_output=True)

ts = np.arange(0, 50, 0.01)
distances = [np.linalg.norm(lorenz_sol2.sol(t) - lorenz_sol1.sol(t)) for t in ts]
distances = np.array(distances)
distances = np.maximum(distances, 1e-16)  # 避免log(0)

plt.figure(figsize=(8, 5))
plt.semilogy(ts, distances, label='距离')
plt.xlabel('$t$')
plt.ylabel('$|\\delta x(t)|$')
plt.title('两条轨迹之间的距离（对数坐标）')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()

我们确实看到在足够大的 $\rho$ 值下，距离在一个显著的窗口内呈指数增长（在半对数图上为直线），然后饱和到吸引子直径的量级。

指数增长的速率被称为李雅普诺夫指数；有更好的方法可以更准确地测量它。

由于洛伦兹方程模拟大气，我们可以预期大气的动力学至少同样复杂。

我们还可以绘制每条轨迹的 $x$ 坐标进行比较：

fig, ax = plt.subplots(figsize=(10, 5))
t_plot = np.linspace(0, 50, 5000)

ax.plot(t_plot, lorenz_sol1.sol(t_plot)[0], label='原始', lw=1)
ax.plot(t_plot, lorenz_sol2.sol(t_plot)[0], label='扰动', lw=1, alpha=0.8)
ax.set_xlabel('$t$')
ax.set_ylabel('$x(t)$')
ax.legend(loc='upper left')
plt.tight_layout()
plt.show()

在足够长的时间后，轨迹分离并表现得截然不同。

现在让我们看看这个模型中的"气候"，即各坐标在长时间内的平均值：

T = 1000.0
rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0

lorenz_sol3 = solve_ivp(lorenz, (0, T), [1 + eps, 1 + eps, 1 + eps],
                        args=(sigma, rho, beta),
                        max_step=0.05, dense_output=True)

lorenz_sol4 = solve_ivp(lorenz, (0, T), [1, 1, 1],
                        args=(sigma, rho, beta),
                        max_step=0.05, dense_output=True)

ts_avg = np.linspace(T/2, T, 10000)

avg3 = np.mean([np.abs(lorenz_sol3.sol(t)) for t in ts_avg], axis=0)
avg4 = np.mean([np.abs(lorenz_sol4.sol(t)) for t in ts_avg], axis=0)

print(f"扰动初始条件的长时间平均绝对值 (x, y, z): {avg3}")
print(f"原始初始条件的长时间平均绝对值 (x, y, z): {avg4}")

扰动初始条件的长时间平均绝对值 (x, y, z): [ 6.45770063  6.91025482 23.54638523]
原始初始条件的长时间平均绝对值 (x, y, z): [ 6.5037865   6.94621232 23.59982518]

我们看到每个分量的平均绝对值大致相同，即使个别轨迹截然不同。这是一个例子，说明即使个体行为非常不同，统计性质也可以相同，这启发了气候——即"平均天气"——可以是稳定的这一观点，即使逐日变化差异很大。

回到开头的问题：如何预测气候？

我们刚才看到，由于混沌（蝴蝶效应），天气系统的长期行为不可预测。那么气候科学家是如何做出预测的呢？

关键在于：气候不是预测具体某天的温度，而是预测温度的统计性质（如长期平均值）。正如我们在洛伦兹系统中看到的，即使个别轨迹截然不同，长期统计量却是稳定的。

现在让我们建立一个最简单的气候模型——零维能量平衡模型。

零维能量平衡模型

地球的温度由能量平衡决定：

$$C \frac{dT}{dt} = \underbrace{\frac{(1-\alpha)S}{4}}_{\text{吸收的太阳辐射}} - \underbrace{(A + BT)}_{\text{向外热辐射}} + \underbrace{a \ln\frac{\mathrm{CO}_2}{\mathrm{CO}_{2,0}}}_{\text{温室效应}}$$

在平衡态（工业革命前温度 $T_0 = 14°C$），吸收的太阳辐射与向外热辐射相等，可消去常数 $A$，得到简化形式：

$$C \frac{dT}{dt} = B(T_0 - T) + a \ln\frac{\mathrm{CO}_2}{\mathrm{CO}_{2,0}}$$

参数	含义	典型值
$C$	热容量（大气+上层海洋）	51 J/(m²·°C)
$B$	气候反馈参数	1.3 W/(m²·°C)
$a$	温室效应系数	5.0 W/m²
$\mathrm{CO}_{2,0}$	工业革命前CO₂浓度	280 ppm
$T_0$	工业革命前全球平均温度	14°C

C = 51.0
B_clim = 1.3
forcing_coef = 5.0
CO2_0 = 280.0
T0 = 14.0

def CO2_model(t):
    return CO2_0 * (1 + (t / 220) ** 3)

def climate_ode(t, T):
    co2 = CO2_model(t)
    greenhouse = forcing_coef * np.log(co2 / CO2_0)
    return (1 / C) * (B_clim * (T0 - T) + greenhouse)

sol_climate = solve_ivp(climate_ode, [0, 170], [T0],
                        t_eval=np.linspace(0, 170, 1000))

plt.figure(figsize=(10, 5))
plt.plot(sol_climate.t + 1850, sol_climate.y[0], 'r-', lw=3)
plt.axhline(y=T0, color='gray', ls='--', alpha=0.5, label=f'工业革命前温度 {T0}°C')
plt.axhline(y=T0 + 2, color='orange', ls='--', alpha=0.5, label=f'巴黎协定目标 +2°C')
plt.xlabel('年份')
plt.ylabel('全球平均温度 (°C)')
plt.title('零维能量平衡模型的温度预测')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

try:
    url = 'https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.csv'
    co2_data = pd.read_csv(url, comment='#')
    has_data = True
except Exception:
    has_data = False

fig, ax1 = plt.subplots(figsize=(10, 5))

t_model = np.linspace(0, 180, 200)
co2_model_vals = [CO2_model(t) for t in t_model]
ax1.plot(t_model + 1850, co2_model_vals, 'r-', lw=2, label='CO2 模型')

if has_data:
    ax1.scatter(co2_data['decimal date'], co2_data['average'],
                s=3, alpha=0.5, color='blue', label='Mauna Loa 实测数据')

ax1.set_xlabel('年份')
ax1.set_ylabel('CO2 浓度 (ppm)')
ax1.set_title('大气 CO2 浓度：模型与实测对比')
ax1.legend()
ax1.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

C_rcp = 100.0
B_rcp = 0.3
alpha_rcp = 1.5
T0_rcp = 14.0

def co2_historical(year):
    return 280 * np.exp(0.0023 * (year - 1850))

def co2_low(year):
    if year <= 2020:
        return co2_historical(year)
    t_f = year - 2020
    return 415 + t_f * 0.5 - t_f**2 * 0.001

def co2_high(year):
    if year <= 2020:
        return co2_historical(year)
    return 415 * np.exp(0.0133 * (year - 2020))

def climate_rcp_ode(t, T, co2_func):
    co2 = co2_func(t)
    forcing = alpha_rcp * np.log(co2 / 280)
    return (1 / C_rcp) * (forcing - B_rcp * (T - T0_rcp))

years = np.linspace(1850, 2100, 500)

sol_low = solve_ivp(lambda t, T: climate_rcp_ode(t, T, co2_low),
                    [1850, 2100], [T0_rcp], t_eval=years)
sol_high = solve_ivp(lambda t, T: climate_rcp_ode(t, T, co2_high),
                     [1850, 2100], [T0_rcp], t_eval=years)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

co2_low_vals = [co2_low(y) for y in years]
co2_high_vals = [co2_high(y) for y in years]

ax1.plot(years, co2_low_vals, '#00aaff', lw=2.5, label='RCP2.6 (低碳)')
ax1.plot(years, co2_high_vals, '#e57347', lw=2.5, label='RCP8.5 (高碳)')
ax1.scatter([2026], [co2_low(2026)], color='black', s=80, zorder=5)
ax1.set_xlabel('年份')
ax1.set_ylabel('CO2 浓度 (ppm)')
ax1.set_title('未来 CO2 排放情景')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2.plot(sol_low.t, sol_low.y[0], '#00aaff', lw=2.5, label='RCP2.6 (低碳)')
ax2.plot(sol_high.t, sol_high.y[0], '#e57347', lw=2.5, label='RCP8.5 (高碳)')
ax2.axhline(y=T0_rcp + 2, color='orange', ls='--', alpha=0.5, label='巴黎协定 +2°C')
ax2.set_xlabel('年份')
ax2.set_ylabel('温度 (°C)')
ax2.set_title('对应温度变化预测')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

小结：天气、气候与混沌

• 天气不可预测：洛伦兹方程展示了确定性混沌——微小的初始差异被指数放大（蝴蝶效应），使得逐日天气预报在一周后失效。
• 气候可以预测：虽然我们不能预测具体某天的温度，但温度的长期统计性质（如全球平均温度）是由能量平衡决定的，可以用相对简单的模型来预测。
• 关键区别：天气是初值问题（对初始条件敏感），气候是边值问题（由能量收支决定）。

我们利用非线性动力学来理解天气预测的局限，同时用能量平衡模型来预测气候变化的趋势。

交互式探索（可选）

    @interact(rho=FloatSlider(min=0, max=100, step=0.1, value=28, description='ρ'))
    def plot_lorenz_interactive(rho=28):
        sigma = 10.0
        beta = 8.0 / 3.0
        
        sol = solve_ivp(lorenz, (0, 50), [0.01, 0.01, 0.01],
                        args=(sigma, rho, beta),
                        max_step=0.02, dense_output=True)
        
        fig = plt.figure(figsize=(8, 6))
        ax = fig.add_subplot(111, projection='3d')
        ax.plot(sol.y[0], sol.y[1], sol.y[2], lw=0.5, alpha=0.8)
        ax.set_xlabel('$x$')
        ax.set_ylabel('$y$')
        ax.set_zlabel('$z$')
        ax.set_xlim(-25, 25)
        ax.set_ylim(-25, 25)
        ax.set_zlim(0, 60)
        ax.set_title(f'洛伦兹吸引子 ($\\rho = {rho:.1f}$)')
        plt.tight_layout()
        plt.show()

interactive(children=(FloatSlider(value=28.0, description='ρ'), Output()), _dom_classes=('widget-interact',))

    @interact(a=FloatSlider(min=0, max=5, step=0.1, value=1.0, description='a'),
              b=FloatSlider(min=0, max=5, step=0.1, value=1.5, description='b'))
    def plot_brusselator_interactive(a=1.0, b=1.5):
        t_span = (0, 10)
        
        fig, ax = plt.subplots(figsize=(7, 7))
        ax.set_facecolor('black')
        
        for x0 in np.arange(0, 5.5, 1.0):
            for y0 in np.arange(0, 5.5, 1.0):
                sol = solve_ivp(brusselator, t_span, [x0, y0], args=(a, b),
                                max_step=0.02, dense_output=True)
                ax.plot(sol.y[0], sol.y[1], lw=1.2, alpha=0.8)
        
        ax.set_xlim(0, 5)
        ax.set_ylim(0, 5)
        ax.set_aspect('equal')
        ax.set_xlabel('$x$')
        ax.set_ylabel('$y$')
        ax.set_title(f'布鲁塞尔子模型 ($a={a:.1f}, b={b:.1f}$)')
        plt.tight_layout()
        plt.show()

interactive(children=(FloatSlider(value=1.0, description='a', max=5.0), FloatSlider(value=1.5, description='b'…

本讲小结

1. ODE回顾：解析解（常数增长、指数趋近、含外力）与数值解（欧拉法、solve_ivp）
2. 一维动力学：不动点及其稳定性、向量场与相图、分岔与双稳态、滞后效应
3. 二维动力学：布鲁塞尔子模型、Hopf分岔与极限环振荡
4. 三维混沌：洛伦兹方程、蝴蝶效应、李雅普诺夫指数
5. 气候模型：零维能量平衡模型、温室效应、未来气候情景预测