{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cf321755-f33f-47b8-9aaf-e3aec2917fda",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# 程序设计与计算思维\n",
    "\n",
    "Computer Programming and Computational Thinking\n",
    "\n",
    "## 第 21 讲：一维的平流与扩散\n",
    "\n",
    "2025—2026学年度春季学期\n",
    "\n",
    "清华大学 韩文弢"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df7a1cea-50fb-4b4e-b64c-d96e1fafd00c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## 问题引入\n",
    "\n",
    "前面研究的是标量在时间维度上的动力学，例如地球温度如何随时间变化。\n",
    "\n",
    "但地球并不具有单一、均匀的温度；更确切地说，在某一特定时刻，地球上不同地方的温度各不相同，这些不同的温度由于多种机制而随时间变化。\n",
    "\n",
    "本讲将研究两种基本机制：**平流**和**扩散**。考虑海洋中的温度。\n",
    "\n",
    "- 由于海洋是运动中的流体，一团温暖的海水可以因水本身的物理运动而流到新的位置；这就是**平流** (advection)。\n",
    "- 即使水不运动，温度或溶解在流体中的某种物质的较高浓度也可以通过分子机制扩散；这就是**扩散** (diffusion)。\n",
    "\n",
    "这些因素的引入将问题从一维（只有时间）扩展到多维（时间和空间）。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5700c12-4eb2-4c92-be20-0b5aa5be0573",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## 准备工作"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbe657f5-98ef-455c-a5a5-0602747c83cd",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact, Checkbox, FloatSlider, IntSlider\n",
    "\n",
    "plt.rcParams['font.family'] = ['Noto Sans', 'Noto Sans CJK SC']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c350737-dd36-4827-959b-4c16b259c380",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## 时间与空间共同演进"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc939e51-f561-4823-9c69-bba06651c725",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "source": [
    "本讲将限制在一维空间，把温度 $T$ 视为以下函数：\n",
    "\n",
    "$$T(t, x)$$\n",
    "\n",
    "依赖于两个独立变量：\n",
    "\n",
    "- 时间，$t$\n",
    "- 空间，$x$\n",
    "\n",
    "需要为每一对可能的值 $(t, x)$ 计算温度 $T$ 的值，即对于所有时间（$>0$）和所有位置。\n",
    "\n",
    "给定点的温度将因不同的物理过程而变化。需要通过写出描述每个物理过程及其对温度影响的方程来建模这一点。由于现在有**两个**独立变量 $t$ 和 $x$，可以预期最终会得到关于这两个变量的导数，因此给定点的温度时间变化率也依赖于温度在空间中的梯度。这将导出一个**偏微分方程** (partial differential equation, PDE)，关联 $T$ 的**偏导数**。\n",
    "\n",
    "在气候建模的背景下，可以把 $x$ 视为**纬度**，假设同一纬度的所有点具有相同温度。这样就可以建模极地寒冷、赤道温暖的事实，以及热量如何从热处流向冷处。\n",
    "\n",
    "然而，显然无法以此方式建模实际的洋流，那将需要两个甚至三个空间维度。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ec94e29-fc22-4729-a734-7240408082eb",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## 温度剖面与离散化\n",
    "\n",
    "常微分方程需要每个变量的初始值。偏微分方程则需要一个**初始函数** $T_0(x)$ 给出每个位置 $x$ 处的温度。假设位置限制在区间 $[0, L_x]$ 内。\n",
    "\n",
    "为了在计算机上表示连续函数 $T_0(x)$，需要以某种方式进行**离散化** (discretization)，即将连续函数近似为计算机中的有限数值集合。\n",
    "\n",
    "最简单的离散化方法是离散**网格点**（或**节点**）$x_i$（$i = 1, \\dots, N_x$）处对函数进行**采样**。为简单起见，取等间距排列，间距 $x_{i+1} - x_i = \\delta x = L_x / N_x$。\n",
    "\n",
    "可以将网格点取在每个区间的**中心**，这样就有 $N_x$ 个区间和 $N_x$ 个网格点，从 $x_1 = \\delta x/2$ 开始，到 $x_N = L_x - \\delta x / 2$ 结束。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adeada81-c37e-4dd3-8174-aa865f1f5988",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def T0_func(x):\n",
    "    return np.sin(2 * np.pi * x) + 2 * np.cos(4 * np.pi * x) + 0.2\n",
    "\n",
    "\n",
    "Nx = 20\n",
    "Lx = 1.0\n",
    "dx = Lx / Nx\n",
    "xs = np.arange(dx / 2, Lx, dx)\n",
    "\n",
    "print(f'网格点数 Nx = {Nx}')\n",
    "print(f'空间步长 δx = {dx}')\n",
    "print(f'网格点 xs = {np.round(xs, 4)}')\n",
    "\n",
    "x_fine = np.linspace(0, Lx, 1000)\n",
    "T0_fine = T0_func(x_fine)\n",
    "T0_sampled = T0_func(xs)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6),\n",
    "                                gridspec_kw={'height_ratios': [4, 1]})\n",
    "\n",
    "ax1.plot(x_fine, T0_fine, lw=2, label='T₀(x)')\n",
    "ax1.scatter(xs, T0_sampled, c='red', zorder=5, label='采样点')\n",
    "ax1.scatter(xs, np.zeros_like(xs), c='black', alpha=0.5, s=15, label='x 节点')\n",
    "\n",
    "for i in range(len(xs)):\n",
    "    ax1.plot([xs[i] - dx/2, xs[i] + dx/2], [T0_sampled[i], T0_sampled[i]],\n",
    "             c='green', lw=4)\n",
    "    ax1.plot([xs[i] - dx/2, xs[i] - dx/2], [0, T0_sampled[i]],\n",
    "             c='green', lw=1, ls='--', alpha=0.3)\n",
    "    ax1.plot([xs[i] + dx/2, xs[i] + dx/2], [0, T0_sampled[i]],\n",
    "             c='green', lw=1, ls='--', alpha=0.3)\n",
    "\n",
    "for x in xs:\n",
    "    ax1.plot([x, x], [0, T0_func(x)], ls='--', c='black', alpha=0.3)\n",
    "\n",
    "ax1.axhline(0, ls='--', c='gray', alpha=0.5)\n",
    "ax1.set_ylabel('T₀(x)')\n",
    "ax1.legend(loc='upper right')\n",
    "\n",
    "ax2.imshow(T0_sampled.reshape(1, -1), aspect='auto', cmap='RdBu_r',\n",
    "           extent=[xs[0]-dx/2, xs[-1]+dx/2, 0, 1],\n",
    "           vmin=-3, vmax=3)\n",
    "ax2.set_yticks([])\n",
    "ax2.set_xlabel('x')\n",
    "ax2.set_ylabel('热图')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf172d00-8bee-4e3c-9b3e-775afe89f8c3",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## 平流\n",
    "\n",
    "现在把这个**剖面** (profile) 视为代表流体中每个小体积或“流体微团”的温度。假设流体以恒定、均匀的速度 $U$ 向右运动（**均匀**意味着流体各部分的速度相同）。那么温度剖面也应该**随流体移动**。这种随流体一起运动的量（如温度）被称为**示踪物** (tracer)。\n",
    "\n",
    "如果将注意力固定在空间中的某个点，比如网格点 $x_i$，那里的温度将随时间变化，因为流体从旁边流过。温度如何随时间变化取决于相邻网格点的值，因为它们决定了有多少热量被传输**进入**和**离开**当前单元。\n",
    "\n",
    "补充知识：将注意力固定在空间某一点的观点称为**欧拉**观点。另一种选择是跟随流体微团在空间中运动；这称为**拉格朗日**观点。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85e772d4-1083-44de-bc2e-c3151f9b7496",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### 可视化流体中的通量\n",
    "\n",
    "随着流体流过网格点（更确切地说，是以网格点为中心的单元），可视化会发生什么。可视化在流体中运动的示踪粒子：\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adf6c7c0-f113-46e6-952c-042ac70309c2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "N_particles = 5000\n",
    "np.random.seed(42)\n",
    "xx_particles = np.abs(-2 + 4 * np.random.rand(N_particles))**2 - 1.5\n",
    "yy_particles = np.random.rand(N_particles)\n",
    "delta_p = 0.8\n",
    "U_p = 0.2\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "ax.fill([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], alpha=0.5, color='lightblue')\n",
    "ax.scatter(xx_particles, yy_particles, s=1.5, alpha=0.1, c='gray')\n",
    "ax.plot([-1.5, 2], [0, 0], c='black')\n",
    "ax.plot([-1.5, 2], [1, 1], c='black')\n",
    "ax.set_xlim(-2, 2)\n",
    "ax.set_ylim(-0.1, 1.1)\n",
    "ax.set_xlabel('x')\n",
    "ax.set_title('流体中的示踪粒子（t = 0）')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98e0fb01-36ee-4226-b7a1-ab64cc040ffa",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_particles(t=0.0, show_particles=False):\n",
    "    fig, ax = plt.subplots(figsize=(10, 3))\n",
    "\n",
    "    ax.fill([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], alpha=0.5, color='lightblue')\n",
    "\n",
    "    new_xx = xx_particles + U_p * t\n",
    "    ax.scatter(new_xx, yy_particles, s=1.5, alpha=0.1, c='gray')\n",
    "\n",
    "    if show_particles:\n",
    "        cs = np.zeros(N_particles, dtype=int)\n",
    "        for i in range(N_particles):\n",
    "            x = xx_particles[i]\n",
    "            if -U_p * delta_p < x < 0:\n",
    "                cs[i] = 1\n",
    "            elif 1 - U_p * delta_p < x < 1:\n",
    "                cs[i] = 2\n",
    "\n",
    "        mask = cs != 0\n",
    "        ax.scatter(new_xx[mask], yy_particles[mask], s=1.5, alpha=0.5, c=cs[mask])\n",
    "\n",
    "    ax.plot([-1.5, 2], [0, 0], c='black')\n",
    "    ax.plot([-1.5, 2], [1, 1], c='black')\n",
    "    ax.set_xlim(-2, 2)\n",
    "    ax.set_ylim(-0.1, 1.1)\n",
    "    ax.set_xlabel('x')\n",
    "    ax.set_title(f'流体中的示踪粒子（t = {t:.2f}）')\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "interact(plot_particles,\n",
    "         t=FloatSlider(min=0, max=2, step=0.01, value=0, description='t'),\n",
    "         show_particles=Checkbox(value=False, description='显示进入/离开的粒子'));\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b168e1a-57aa-4cdb-bff0-a7f13808c5e0",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### 时间步进\n",
    "\n",
    "希望建模温度剖面如何因流体流动而随时间变化。方法是观察每个单元，询问在给定时间步长 $\\delta t$ 内有多少热量进入和离开该单元。\n",
    "\n",
    "将 $T^n_i$ 称为位置 $x_i$ 处、第 $n$ 个时间步 $t_n$ 的近似（未知）平均值，即 $T(t_n, x_i)$ 的近似，其中 $t_n = n \\, \\delta t$。\n",
    "\n",
    "于是 $T^{n+1}_i \\simeq T(t_n + \\delta t, x_i)$，$T^{n}_{i+1} \\simeq T(t_n, x_i + \\delta x)$。\n",
    "\n",
    "注意：这些算法中的上标 $n$ 并不表示幂；它只是时间步的标签。\n",
    "\n",
    "假设流体以速度 $U$ 向右运动。在持续 $\\delta t$ 的时间步内，单元 $i$ 的温度 $T^n_i$ 因两个原因而改变：\n",
    "\n",
    "- 一些热量**进入**单元 $i$\n",
    "- 一些热量**离开**单元 $i$\n",
    "\n",
    "要计算有多少热量进入和离开，注意只有在距单元边界距离 $U \\, \\delta t$ 以内的流体区域中的热量才会穿越边界。因此，单元 $i$ 中热量的比例为 $(U \\, \\delta t) / \\delta x$ 会穿越边界。\n",
    "\n",
    "因此，有 $T^n_i (U \\delta t) / \\delta x$ 的量将离开单元 $i$ 并进入单元 $i+1$（右侧单元）。类似地，$T^n_{i-1} (U \\delta t) / \\delta x$ 的量将从左侧相邻单元 $i-1$ 进入单元 $i$。\n",
    "\n",
    "由此得到：\n",
    "\n",
    "$$T^{n+1}_i = T^{n}_i + (T^n_{i-1} - T^n_{i})\\, U \\, \\delta t / \\delta x.$$\n",
    "\n",
    "注意：等式右边是时间步 $n$ 的量，左边是时间步 $n+1$ 的量。于是得到了如何从第 $n$ 步**更新**到第 $n+1$ 步的方法。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2467d562-0ac8-4dae-9598-f0f81f895ec2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### 连续极限：平流方程 PDE\n",
    "\n",
    "重新排列前一个方程得到\n",
    "\n",
    "$$\\frac{T^{n+1}_i - T^{n}_i}{\\delta t} = \\frac{T^n_{i-1} - T^n_{i}}{\\delta x}\\, U.$$\n",
    "\n",
    "当 $\\delta t \\to 0$ 和 $\\delta x \\to 0$ 时取连续极限，可以看到是**偏导数**的定义。用 $\\partial$ 表示这些偏导数，得到**平流方程**：\n",
    "\n",
    "$$\\frac{\\partial T(t, x)}{\\partial t} = -U \\frac{\\partial T(t, x)}{\\partial x},$$\n",
    "\n",
    "简写为\n",
    "\n",
    "$$\\frac{\\partial T}{\\partial t} = -U \\frac{\\partial T}{\\partial x}.$$\n",
    "\n",
    "由于 $T$ 是 $x$ 和 $t$ 的函数，且此方程涉及关于两个独立变量的偏导数，因此这是一个**偏微分方程**（PDE）。它描述了函数 $T(t, x)$ 如何随时间和空间连续变化。\n",
    "\n",
    "虽然有一些解析方法可以求解 PDE，但通常需要使用数值方法。这里将研究求解此类方程的简单数值方法。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c939325f-d66e-4738-ba01-c8688b25133e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### 平流方程的数值方法\n",
    "\n",
    "回到以下形式的方程，其中下一时间步的值已被分离：\n",
    "\n",
    "$$T^{n+1}_i = T^{n}_i - \\left( U \\frac{\\delta t}{\\delta x} \\right) (T^n_{i} - T^n_{i-1}).$$\n",
    "\n",
    "在右边最后一项中，需要同一时间步中不同位置 $T$ 值的某种组合。\n",
    "\n",
    "最简单的方法是直接转录向量的第 $i$ 个分量的方程。将 `T` 称为当前向量，即 $\\mathbf{T}^n := (T^n_i)_{i=1, \\ldots, N_x}$，`T_new` 为下一个时间步的新向量：\n",
    "\n",
    "```\n",
    "T_new[i] = T[i] + δt * U * (T[i-1] - T[i]) / δx\n",
    "```\n",
    "\n",
    "但这时出现一个问题：当 $i=0$ 时怎么办？这将试图访问向量 `T` 的索引 $-1$，该索引不存在。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec44a066-472e-4f03-a3f9-cb95b36af972",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### 边界条件\n",
    "\n",
    "这说明必须选择**边界条件**来指定区域边缘发生什么。\n",
    "\n",
    "为简单起见，选择使用**周期性边界条件**。这是一种方便的数学处理，通过将系统环绕在一个环面上，使所有单元处于同等地位，从而使单元 $i=0$ 和 $i=N_x-1$ 成为相邻单元。\n",
    "\n",
    "于是可以写成如下形式，单独写 $i=0$ 的情况（也可以考虑取模）。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3edc7e17-c666-46e6-acd4-962978d8ed93",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "dt = 0.001\n",
    "U_adv = 0.2\n",
    "\n",
    "\n",
    "def advection(T, dt, dx, U):\n",
    "    N = len(T)\n",
    "    T_new = np.zeros_like(T)\n",
    "\n",
    "    for i in range(1, N):\n",
    "        T_new[i] = T[i] - dt * U * (T[i] - T[i-1]) / dx\n",
    "\n",
    "    T_new[0] = T[0] - dt * U * (T[0] - T[N-1]) / dx\n",
    "\n",
    "    return T_new"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a95b0c2-4ee2-49dc-9b52-c25338bb8dd5",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "source": [
    "这是平流方程的单次时间步；它接受当前的 $T$ 向量并返回更新后的新 $T$ 向量。\n",
    "\n",
    "注意，这就像求解 ODE 的欧拉法的一步，但同时更新许多空间坐标。本质上在求解一个耦合的 ODE 组。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77c15afe-655d-46fa-bb47-9c02eee7153f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def temperature_heatmap(x, T, ax=None, vmin=-3, vmax=3):\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(figsize=(10, 1))\n",
    "    ax.imshow(np.array(T).reshape(1, -1), aspect='auto', cmap='RdBu_r',\n",
    "              extent=[x[0]-dx/2, x[-1]+dx/2, 0, 1],\n",
    "              vmin=vmin, vmax=vmax)\n",
    "    ax.set_yticks([])\n",
    "    return ax\n",
    "\n",
    "\n",
    "def evolve(method, xs, dt, params, t_final=10.0, f0=T0_func):\n",
    "    T = np.array([f0(x) for x in xs])\n",
    "\n",
    "    t = 0.0\n",
    "    ts = [t]\n",
    "    results = [T.copy()]\n",
    "\n",
    "    while t < t_final - dt/2:\n",
    "        T_new = method(T, dt, dx, params)\n",
    "        results.append(T_new.copy())\n",
    "        t += dt\n",
    "        ts.append(t)\n",
    "        T = T_new.copy()\n",
    "\n",
    "    return np.array(ts), results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98af02be-1270-4fa2-b2af-161c46a096a8",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "ts_adv, evolution_adv = evolve(advection, xs, dt, U_adv, t_final=2.0)\n",
    "\n",
    "n_snapshot = len(evolution_adv) // 2\n",
    "\n",
    "def plot_advection(n=n_snapshot):\n",
    "    n = min(n, len(evolution_adv) - 1)\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5),\n",
    "                                    gridspec_kw={'height_ratios': [4, 1]})\n",
    "\n",
    "    ax1.scatter(xs, evolution_adv[n], c='blue', zorder=5)\n",
    "    ax1.plot(xs, evolution_adv[n], c='blue', alpha=0.5)\n",
    "    ax1.set_xlim(0, 1)\n",
    "    ax1.set_ylim(-3.1, 3.1)\n",
    "    ax1.set_ylabel('T')\n",
    "    ax1.set_title(f'平流（上游风格）：t = {ts_adv[n]:.2f}')\n",
    "\n",
    "    temperature_heatmap(xs, evolution_adv[n], ax2)\n",
    "    ax2.set_xlabel('x')\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "plot_advection()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bc44bc5-d088-40b9-a60a-90fd3cc20ba4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "interact(plot_advection,\n",
    "         n=IntSlider(min=1, max=len(evolution_adv)-1, step=1, value=1, description='时间步'));\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c6e1784-fdab-40bc-8ac9-3222c802fce8",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "遗憾的是，这种方法并**不能**按预期工作：随着时间的推移，剖面形状没有保持不变，而是在衰减。这是由于近似方式造成的。\n",
    "\n",
    "一种更好的空间导数离散化方法是使用以下**中心差分**：\n",
    "\n",
    "$$\\frac{\\partial T(t_n, x_i)}{\\partial x} \\simeq \\frac{T^n_{i+1} - T^n_{i-1}}{2 \\delta x}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72c5997f-3c82-4ca7-aa17-837b2dc4b1ec",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def advection2(T, dt, dx, U):\n",
    "    N = len(T)\n",
    "    T_new = np.zeros_like(T)\n",
    "\n",
    "    for i in range(1, N-1):\n",
    "        T_new[i] = T[i] - dt * U * (T[i+1] - T[i-1]) / (2 * dx)\n",
    "\n",
    "    T_new[0] = T[0] - dt * U * (T[1] - T[N-1]) / (2 * dx)\n",
    "    T_new[N-1] = T[N-1] - dt * U * (T[0] - T[N-2]) / (2 * dx)\n",
    "\n",
    "    return T_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f25487b-dafc-4951-9529-f5a284c68e48",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "ts_adv2, evolution_adv2 = evolve(advection2, xs, dt, 0.1, t_final=2.0)\n",
    "\n",
    "n_snapshot2 = len(evolution_adv2) // 2\n",
    "\n",
    "def plot_advection2(n=n_snapshot2):\n",
    "    n = min(n, len(evolution_adv2) - 1)\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5),\n",
    "                                    gridspec_kw={'height_ratios': [4, 1]})\n",
    "\n",
    "    ax1.scatter(xs, evolution_adv2[n], c='blue', zorder=5)\n",
    "    ax1.plot(xs, evolution_adv2[n], c='blue', alpha=0.5)\n",
    "    ax1.set_xlim(0, 1)\n",
    "    ax1.set_ylim(-3.1, 3.1)\n",
    "    ax1.set_ylabel('T')\n",
    "    ax1.set_title(f'平流（中心差分）：t = {ts_adv2[n]:.2f}')\n",
    "\n",
    "    temperature_heatmap(xs, evolution_adv2[n], ax2)\n",
    "    ax2.set_xlabel('x')\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "plot_advection2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9c30def-881e-4ed3-9894-eb1bbad4a9d8",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "interact(plot_advection2,\n",
    "         n=IntSlider(min=1, max=len(evolution_adv2)-1, step=1, value=1, description='时间步'));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4abf824f-3eec-4fe9-bb36-1580cc789955",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## 扩散\n",
    "\n",
    "另一个关键的物理过程是**扩散**，可以描述温度或质量如何从高温或高浓度区域向低温或低浓度区域扩散。\n",
    "\n",
    "### 物理机制：随机游走\n",
    "\n",
    "扩散背后的物理机制是**随机运动**：这是描述一群随机游走者在空间和时间中概率分布演化的方程的连续极限。\n",
    "\n",
    "用已有的符号表示随机游走者群的概率分布满足以下时间演化：\n",
    "\n",
    "$$p^{n+1}_i = \\frac{1}{2}(p^n_{i-1} + p^n_{i+1})$$\n",
    "\n",
    "如果游走者只以一定概率跳跃，大概率停留在原地，且这些随机游走者是热量的载体，那么就可以得到\n",
    "\n",
    "$$T^{n+1}_i = \\kappa (T^n_{i-1} - 2 T^n_i + T^n_{i+1}).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cc60d7b-3e2a-4cee-b666-6573e58fc9a7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### 连续极限：热传导方程 PDE\n",
    "\n",
    "引入 $\\delta x$ 作为空间离散化，$\\delta t$ 作为时间步长，得到\n",
    "\n",
    "$$T^{n+1}_i = \\kappa \\frac{\\delta t}{\\delta x^2} (T^n_{i-1} - 2 T^n_i + T^n_{i+1}).$$\n",
    "\n",
    "连续极限是以下**热传导方程**或**扩散方程**：\n",
    "\n",
    "$$\\frac{\\partial T}{\\partial t} = \\kappa \\frac{\\partial^2 T}{\\partial x^2}.$$\n",
    "\n",
    "其中 $\\kappa$ 是**热扩散率**，描述热量扩散的速度。在质量扩散的背景下，等价的是**扩散系数** $D$。\n",
    "\n",
    "要获得求解此方程的数值方法，再次需要离散化，特别是二阶导数。一种可能的离散化是\n",
    "\n",
    "$$\\frac{\\partial^2 T}{\\partial x^2}(t_n, x_i) \\simeq \\frac{T^n_{i+1} - 2 T^n_i + T^n_{i-1}}{\\delta x^2}.$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07f316f8-6465-4bd4-b98b-415d60c6f115",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def diffusion(T, dt, dx, D):\n",
    "    N = len(T)\n",
    "    T_new = np.zeros_like(T)\n",
    "\n",
    "    for i in range(1, N-1):\n",
    "        T_new[i] = T[i] + dt * D * (T[i+1] - 2*T[i] + T[i-1]) / (dx**2)\n",
    "\n",
    "    T_new[0] = T[0] + dt * D * (T[1] - 2*T[0] + T[N-1]) / (dx**2)\n",
    "    T_new[N-1] = T[N-1] + dt * D * (T[0] - 2*T[N-1] + T[N-2]) / (dx**2)\n",
    "\n",
    "    return T_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d9920bd-331e-4807-b049-363e60cdc2f7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "ts_diff, evolution_diff = evolve(diffusion, xs, dt, 0.01, t_final=5.0)\n",
    "\n",
    "n_snapshot_diff = len(evolution_diff) // 2\n",
    "\n",
    "def plot_diffusion(n=n_snapshot_diff):\n",
    "    n = min(n, len(evolution_diff) - 1)\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5),\n",
    "                                    gridspec_kw={'height_ratios': [4, 1]})\n",
    "\n",
    "    ax1.scatter(xs, evolution_diff[n], c='blue', zorder=5)\n",
    "    ax1.plot(xs, evolution_diff[n], c='blue', alpha=0.5)\n",
    "    ax1.set_xlim(0, 1)\n",
    "    ax1.set_ylim(-3.1, 3.1)\n",
    "    ax1.set_ylabel('T')\n",
    "    ax1.set_title(f'扩散：t = {ts_diff[n]:.2f}')\n",
    "\n",
    "    temperature_heatmap(xs, evolution_diff[n], ax2)\n",
    "    ax2.set_xlabel('x')\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "plot_diffusion()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7386ae35-ab57-47cd-8aae-c7d8ac7e1823",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "interact(plot_diffusion,\n",
    "         n=IntSlider(min=1, max=len(evolution_diff)-1, step=1, value=1, description='时间步'));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d465da1e-b1ee-4468-bd2b-aca7d6406d6b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## 平流-扩散 PDE\n",
    "\n",
    "最后，可以将两种机制结合起来，描述一个同时被平流（以恒定速度运动）和扩散的示踪物。这基本上利用了平流和扩散函数的复合：\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bf3f9bc-1533-4061-81ad-6e22670bc062",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def advection_diffusion(T, dt, dx, params):\n",
    "    U, D = params\n",
    "    temp = advection2(T, dt, dx, U)\n",
    "    return diffusion(temp, dt, dx, D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55259d24-6e43-4f5e-b524-ed191ecef8ab",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "UU_default = 0.5\n",
    "DD_default = 0.01\n",
    "\n",
    "ts_ad, results_ad = evolve(advection_diffusion, xs, dt,\n",
    "                           (UU_default, DD_default), t_final=2.0)\n",
    "\n",
    "n_snapshot_ad = len(results_ad) // 2\n",
    "\n",
    "def plot_advection_diffusion(U=UU_default, D=DD_default, n=n_snapshot_ad):\n",
    "    ts_local, results_local = evolve(advection_diffusion, xs, dt,\n",
    "                                     (U, D), t_final=2.0)\n",
    "    n = min(n, len(results_local) - 1)\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 5),\n",
    "                                    gridspec_kw={'height_ratios': [4, 1]})\n",
    "\n",
    "    ax1.scatter(xs, results_local[n], c='blue', zorder=5)\n",
    "    ax1.plot(xs, results_local[n], c='blue', alpha=0.5)\n",
    "    ax1.set_xlim(0, 1)\n",
    "    ax1.set_ylim(-3.1, 3.1)\n",
    "    ax1.set_ylabel('T')\n",
    "    ax1.set_title(f'平流-扩散（U={U:.2f}, D={D:.3f}）：t = {ts_local[n]:.2f}')\n",
    "\n",
    "    temperature_heatmap(xs, results_local[n], ax2)\n",
    "    ax2.set_xlabel('x')\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "plot_advection_diffusion()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1a0cd30-668a-498b-af52-de4c867a12e5",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "interact(lambda n: plot_advection_diffusion(n=n),\n",
    "         n=IntSlider(min=1, max=2000, step=1, value=1, description='时间步'));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bec5c35-ee96-4949-acd9-7651e4d4285d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "-"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "interact(plot_advection_diffusion,\n",
    "         U=FloatSlider(min=-1, max=1, step=0.01, value=0.5, description='U'),\n",
    "         D=FloatSlider(min=-0.2, max=0.2, step=0.001, value=0.01, description='D'),\n",
    "         n=IntSlider(min=1, max=2000, step=1, value=1, description='时间步'));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35a8cce1-016e-4217-a848-8d0da464af01",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## 本讲小结\n",
    "\n",
    "本讲从 ODE（常微分方程）迈向 PDE（偏微分方程），引入了空间维度的离散化与建模。从**计算思维**的角度：\n",
    "\n",
    "**1. 抽象与维度的扩展**\n",
    "\n",
    "此前，温度 $T$ 仅是时间 $t$ 的函数（ODE）。本讲将其扩展为 $T(t, x)$，依赖时间和空间两个独立变量，由此引入**偏微分方程**。这是计算思维中\"抽象层次提升\"的典型例子：同一个时间步进框架，从处理标量变为处理向量（空间网格上的离散值），代码结构几乎不变。\n",
    "\n",
    "**2. 离散化（Discretization）**\n",
    "\n",
    "连续函数 $T(x)$ 无法直接在计算机中表示。本讲采用**网格点采样**的方法，将空间区间 $[0, L_x]$ 均分为 $N_x$ 个单元，在每个单元中心存储一个温度值。这种空间离散化与第17讲的时间离散化（欧拉法）相结合，形成了求解 PDE 的基础：\n",
    "\n",
    "- **空间离散化**：上游差分 vs 中心差分，以及二阶导数的离散化\n",
    "- **时间离散化**：前向欧拉法从 ODE 推广到 PDE\n",
    "\n",
    "**3. 分解与组合（Decomposition & Composition）**\n",
    "\n",
    "平流和扩散是两个独立的物理过程，各自对应一个 PDE：\n",
    "- 平流方程：$\\partial T / \\partial t = -U \\, \\partial T / \\partial x$\n",
    "- 扩散方程：$\\partial T / \\partial t = \\kappa \\, \\partial^2 T / \\partial x^2$\n",
    "\n",
    "代码中分别实现了 `advection2()` 和 `diffusion()` 函数，然后通过**函数复合** `advection_diffusion()` 将两者组合。这种\"先分解、再组合\"的策略使得代码模块化，且易于扩展新的物理过程。\n",
    "\n",
    "**4. 建模：从物理直觉到代码**\n",
    "\n",
    "- 平流：流体携带着温度移动——对应**上游差分**格式\n",
    "- 扩散：连接到第13讲的随机游走——随机游走的连续极限就是扩散方程\n",
    "\n",
    "这两个物理过程的数值实现都直接来源于对物理机制的理解，体现了\"物理 → 数学 → 代码\"的建模路径。\n",
    "\n",
    "扩展：根据本讲内容建立地球的一维能量平衡模型，计算得出赤道热、两极冷的温度分布。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.3"
  },
  "rise": {
   "scroll": true,
   "slideNumber": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}