MIT-18-S191-计算机思维导论笔记-全-

MIT 18.S191 计算机思维导论笔记（全）

L2：图像处理与变换 🖼️

在本节课中，我们将学习如何对图像进行各种变换。我们将从简单的像素操作开始，逐步深入到线性组合和卷积等核心数学概念。通过这些操作，我们不仅能处理图像，还能理解背后广泛应用的数学和计算机科学思想。

课程概述

本节课将介绍图像处理的基本操作，包括像素采样、线性组合以及卷积。我们将使用Julia语言进行演示，并解释这些操作背后的数学原理。

图像像素化与上采样

首先，我们来看看如何通过下采样使图像像素化，以及如何通过上采样恢复图像。

下采样是通过每隔r个像素取一个点来实现的。当r值增大时，图像会变得更加像素化，因为获取的像素点变少了。当r为1时，我们得到原始图像。

上采样则通过一种称为克罗内克积的数学操作实现。这个操作将每个像素替换为一个r乘r的像素块，块中所有像素的颜色都与原像素相同。这样就能快速放大图像。

线性组合与图像混合

上一节我们介绍了像素操作，本节中我们来看看如何将图像进行线性组合。线性组合是数学、应用数学和工程学中无处不在的概念。

线性组合的核心是两个基本操作：缩放对象和组合多个对象。在图像处理中，这意味着我们可以调整图像的亮度，然后将多张图像叠加在一起。

以下是缩放图像亮度的操作：

# 将图像亮度乘以常数c
brightened_image = c .* original_image

当c大于1时，图像变亮；当c小于1时，图像变暗。但需要注意的是，颜色值有上限，过度增亮会导致颜色饱和变为白色。

接下来，我们可以组合多张图像。例如，将一张正常图像和一张上下颠倒的图像进行混合。

以下是创建图像线性组合的示例：

# alpha控制混合比例
combined_image = alpha .* rightsideup_image + (1-alpha) .* upsidedown_image

当alpha为1时，只显示正常图像；为0时，只显示颠倒图像；为0.5时，两者等比例混合。有趣的是，人脑对正立面孔的识别优于倒立面孔，这是一种被称为“面孔倒置效应”的知觉现象。

如果线性组合中所有权重均为非负且总和为1，这种组合被称为凸组合。

卷积：图像处理的核心

上一节我们探索了线性组合，本节中我们来看看图像处理中更强大的工具：卷积。卷积在机器学习等领域有广泛应用，是计算机识别图像的关键。

卷积操作是将一个图像与一个称为“核”的小矩阵进行运算。核在图像上滑动，每个位置的输出是核覆盖区域内像素值的加权和。

以下是几种常见的卷积核及其效果：

• 模糊核：取像素及其周围像素的平均值，使图像变模糊。
• 边缘检测核：突出显示图像中颜色或亮度急剧变化的区域，即边缘。
• 锐化核：通常是恒等核和边缘检测核的线性组合，在保留原图的同时增强边缘。
• 梯度核：用于近似计算图像在x或y方向上的导数（即变化率）。

卷积的计算复杂度大致为图像像素数乘以核的大小。虽然小核计算更快，但现代GPU（图形处理器）因其并行架构，可以极其高效地处理这种规律的卷积运算，这也是它们被广泛用于机器学习和图像处理的原因。

高斯模糊

在众多模糊方法中，高斯模糊因其良好的数学特性而备受青睐。与简单的平均模糊不同，高斯模糊核的权重不是均匀的，而是中心高、四周低，符合二维高斯分布。

高斯核的权重由以下公式决定：

# 二维高斯函数近似
weight = exp(-(x^2 + y^2) / 2)

以下是高斯模糊的特点：

• 它更重视中心像素，对远处像素影响较小，因此能产生更平滑、更自然的模糊效果。
• 核越大，模糊程度越高。
• 在软件中，可以通过定义偏移数组来灵活设定核的中心位置，使其在处理边界时更加方便。

离散与连续

在图像处理中，我们处理的是离散的像素点。然而，许多操作在连续数学中都有对应的概念。

例如，之前提到的x方向梯度核，本质上就是离散版本的导数。这提醒我们，离散数学和连续数学之间的鸿沟并不像有时看起来那么大。许多核心概念是相通的，只是在不同领域有不同的表现形式和名称。在学习中，留意这些联系会加深我们的理解。

边界处理

在进行卷积时，当核移动到图像边缘时，会出现“越界”问题。程序员需要决定如何处理这些边界情况。

常见的边界处理策略包括：

• 假设边界外的像素值为0（补零）。
• 复制边缘的像素值。
• 镜像图像边缘的像素。
• 环绕图像（即认为图像是循环的）。

不同的应用场景可能需要不同的边界处理方式，优秀的软件会提供灵活的选择。

课程总结

本节课中我们一起学习了图像处理的基本变换。我们从像素采样和上/下采样开始，理解了如何改变图像分辨率。接着，我们探讨了线性组合，用它来混合和调整图像，并接触了凸组合的概念。然后，我们深入学习了卷积这一强大工具，用它实现了模糊、锐化、边缘检测等效果，并了解了高斯模糊的原理。最后，我们讨论了离散与连续数学的联系，以及卷积中边界处理的重要性。通过这些内容，我们不仅学会了处理图像，更看到了背后广泛应用的数学思想。

L3：变换与自动微分 🧮

在本节课中，我们将要学习变换（Transformations）与自动微分（Automatic Differentiation）的核心概念。我们将从简单的函数定义开始，逐步深入到向量函数、线性变换，并探索计算机如何高效地计算导数。通过直观的图像变换演示和代码示例，你将建立起对线性代数与多元微积分的新直觉。

函数定义：三种方式 📝

上一节我们介绍了课程的整体目标，本节中我们来看看在Julia中定义函数的几种不同方式。每种方式都有其适用场景，理解它们有助于编写更清晰、更高效的代码。

以下是三种主要的函数定义方式：

• 短格式：适用于简单的单行函数。例如，定义平方函数：

  f₁(x) = x^2
  

• 匿名格式：函数本身没有名称，通常用于作为参数传递。例如，定义正弦函数：

  x -> sin(x)
  

  你可以将其赋值给一个变量来使用：

  my_sin = x -> sin(x)
  

• 长格式：使用 function 和 end 关键字，适合多行或逻辑复杂的函数。可以设置默认参数。

  function f₃(x; α=3)
      return x^α
  end
  

自动微分入门 ⚡

了解了函数的基本定义后，我们进入本节课的核心技术之一：自动微分。这与你在微积分课上学到的符号求导或数值近似求导都不同。

我们使用 ForwardDiff 包来进行自动微分。

• 标量函数求导：对于单变量函数，求导变得非常简单。

  using ForwardDiff
  f₁(x) = x^2
  derivative = ForwardDiff.derivative(f₁, 5) # 在x=5处的导数
  

• 梯度计算：对于多变量标量函数，自动微分可以计算其梯度（即所有偏导数组成的向量）。

  f₅(x, y, z) = 5 * sin(x*y) + 2*y / (4*z)
  # 在点(1,2,3)计算梯度
  grad = ForwardDiff.gradient(v -> f₅(v[1], v[2], v[3]), [1, 2, 3])
  

  从数学上看，梯度 ∇f 给出了函数增长最快的方向。

向量函数与变换 🌀

上一节我们处理了输入为向量、输出为标量的函数，本节中我们来看看输入和输出都是向量的函数，这通常被称为变换。

以下是几个变换的例子，包括线性和非线性的：

• 恒等变换：最简单的线性变换。

  identity_transform((x, y)) = [x, y]
  

• 缩放变换：沿x轴或y轴进行缩放。

  scale_x(α) = ((x, y),) -> [α*x, y]
  scale_y(α) = ((x, y),) -> [x, α*y]
  

• 旋转变换：将点绕原点旋转角度θ。

  rotate(θ) = ((x, y),) -> [cos(θ)*x - sin(θ)*y, sin(θ)*x + cos(θ)*y]
  

• 剪切变换：

  shear(α) = ((x, y),) -> [x + α*y, y]
  

• 非线性扭曲变换：旋转角度取决于点到原点的距离。

  warp(α) = ((x, y),) -> rotate(α * sqrt(x^2 + y^2))((x, y))
  

雅可比矩阵：变换的导数 📐

对于向量到向量的变换，其“导数”是一个矩阵，称为雅可比矩阵。它描述了变换在局部的最佳线性近似。

对于变换 F([x, y]) = [u(x,y), v(x,y)]，其雅可比矩阵 J 为：

J = [ ∂u/∂x  ∂u/∂y ]
    [ ∂v/∂x  ∂v/∂y ]

我们可以用自动微分轻松计算雅可比矩阵：

F(p) = warp(3)(p) # 以α=3的扭曲变换为例
J = ForwardDiff.jacobian(F, [4, 5])

雅可比矩阵的行列式绝对值给出了该点附近面积的局部缩放因子。

矩阵即线性变换 🖼️

现在，让我们将理论可视化。一个矩阵本质上代表了一个线性变换。矩阵乘法 A * v 就是将线性变换 A 应用于向量 v。

所有线性变换都能保持网格线的“平行四边形”结构。如果一个变换能将均匀的方格网格映射为另一组均匀的平行四边形网格，那么它就可以用一个矩阵来表示。

• 行列式的几何意义：矩阵 A 的行列式 det(A) 表示该线性变换导致的面积缩放比例。如果 det(A) = 1，则变换是保面积的（如旋转、剪切）。

对于非线性变换，虽然整体网格会扭曲，但在任意点的无穷小邻域内，它仍然由一个雅可比矩阵（即一个线性变换）来近似描述。

总结 🎯

本节课中我们一起学习了：

1. Julia中定义函数的多种方式。
2. 自动微分的原理与应用，它不同于符号或数值求导，是机器学习等领域的关键技术。
3. 向量值函数，即变换，包括线性和非线性变换的例子。
4. 雅可比矩阵作为多变量变换的导数，描述了局部线性行为。
5. 将矩阵理解为线性变换，并通过图像变换直观感受行列式等概念的几何意义。

通过结合编程实践与数学直观，我们希望你能对线性代数和多元微积分产生新的、更深刻的理解。

L4：变换：线性变化与非线性变化，组合变换 🌀

在本节课中，我们将要学习图像变换的核心概念，包括线性变换与非线性变换，以及如何组合这些变换。我们将通过直观的图像操作来理解抽象的数学概念，并探索函数与软件的组合性。

---

概述

上一节我们介绍了Julia编程的基本操作和图像处理作为学习背景。本节中，我们将深入探讨变换的数学原理。我们将看到，变换不仅仅是移动像素，它背后对应着深刻的线性代数概念。通过动手操作，我们可以直观地理解矩阵如何作用于图像，以及线性与非线性变换的区别。

---

线性变换：从图像到矩阵

线性变换是一类特殊的变换，它满足两个核心性质：缩放性和可加性。在图像上，一个直观的表现是：一个由等大小矩形组成的网格，在经过线性变换后，会变成一个由全等平行四边形组成的网格。

以下是线性变换的数学定义。一个变换 T 是线性的，当且仅当对于任意向量 v₁, v₂ 和任意标量 c₁, c₂，都满足：

T(c₁ * v₁ + c₂ * v₂) = c₁ * T(v₁) + c₂ * T(v₂)

在代码中，我们可以用矩阵乘法来实现一个通用的线性变换：

linear_transform(M) = v -> M * v

其中 M 是一个2x2矩阵，v 是坐标向量 [x, y]。

以下是几种常见的线性变换示例：

• 缩放：改变图像在x或y方向的大小。
• 旋转：围绕原点旋转图像。
• 剪切：使图像在某一方向上发生倾斜，而另一方向的坐标保持不变。
• 翻转：沿x轴或y轴镜像图像。

所有这些变换都可以通过一个2x2矩阵来定义。矩阵的每一列，实际上代表了变换后标准基向量 [1, 0] 和 [0, 1] 的新位置。

---

非线性变换：打破网格的规则

并非所有变换都是线性的。非线性变换不满足上述的线性条件，它们在图像上会产生更复杂、更“弯曲”的效果。一个简单的例子是平移（移动整个图像），因为它不满足“零向量映射到零向量”这一线性变换的基本要求。

在非线性变换下，矩形网格可能会变成各种曲线网格。然而，有一个重要的观察：在非常小的局部区域内，任何光滑的非线性变换都近似于一个线性变换。这是许多应用数学方法的基础。

以下是一些非线性变换的例子：

• 平移：(x, y) -> (x + a, y + b)
• 非线性剪切：例如 (x, y) -> (x + α*y^2, y)
• 极坐标变换：将图像从直角坐标系映射到极坐标系。
• 自定义扭曲函数：例如 (x, y) -> (x + α*sin(5x), y)

你可以尝试输入任何数学函数来创建自己的变换，观察图像如何被扭曲，这是探索函数行为的绝佳方式。

---

透视变换：一个有趣的边界案例

一个常见的误解是，所有“将直线映射为直线”的变换都是线性的。透视变换（如模拟人眼观看三维场景的投影）就是一个反例。在透视变换中，原本平行的直线在图像中会相交于“消失点”，且等大的矩形会变成大小不一的梯形，这违反了线性变换中“全等平行四边形”的特性。

因此，“保持直线性”不足以定义线性变换。透视变换在数学上通常用齐次坐标和3x3矩阵（对于二维）或4x4矩阵（对于三维）来表示，并在最后一步进行除法，这使得它整体上是非线性的。

---

变换的组合与逆变换

组合是数学和计算机科学中的一个核心思想。组合变换意味着按顺序应用多个变换。在Julia中，我们可以使用函数复合运算符 ∘。

例如，先旋转再剪切可以表示为：

composed_transform = shear(α) ∘ rotate(θ)

一个特别重要的概念是逆变换。如果变换 T 将点A移动到点B，那么它的逆变换 T⁻¹ 恰好将点B移回点A。当组合一个变换和它的逆变换时，结果将是恒等变换（什么都不做）：

T⁻¹ ∘ T = Identity

在图像处理中，这意味着先应用一个扭曲，再应用其反向扭曲，图像将恢复原状。

---

代码风格与函数设计

在编写变换函数时，代码的可读性很重要。比较以下两种定义缩放函数的方式：

# 方式一：使用向量索引
scale_1(α) = v -> [α*v[1], v[2]]

# 方式二：解构参数，命名更清晰
scale_2(α) = ((x, y),) -> (α*x, y)

第二种方式直接使用 x 和 y 作为参数名，更清晰地表达了函数的意图。

当变换需要参数时（如旋转角度α），我们使用匿名函数格式来定义，这允许我们轻松创建一系列参数化的变换函数库。

---

软件的组合性

组合性的思想不仅限于数学函数，它也是优秀软件设计的目标。在理想情况下，不同人编写的代码模块应该能够无缝地协同工作，而无需预先进行大量协调。例如：

• 一个用于求解微分方程的模块。
• 一个用于优化参数的模块。
• 一个用于机器学习的神经网络模块。

如果这些模块设计良好（具有清晰的接口和一致的行为），我们就可以像组合数学函数一样将它们组合起来，用优化器包裹求解器，或用神经网络作为微分方程的一部分。Julia语言社区正致力于实现这种高度的可组合性，让科研代码像乐高积木一样易于拼接。

---

总结

本节课中我们一起学习了变换的核心概念。我们理解了线性变换可以通过矩阵乘法实现，并将网格变为平行四边形；而非线性变换则能产生更丰富、更复杂的扭曲效果。我们探讨了透视变换作为非线性变换的例子，并强调了“直线映射到直线”并非线性的充分条件。最后，我们深入了解了变换的组合与逆变换的概念，并将这种组合性提升到了软件设计哲学的高度。

通过动手操作图像，我们将抽象的代数定义转化为直观的视觉体验，这正是计算思维的魅力所在：利用计算机的交互能力，去探索和深化我们对数学世界的理解。

L5：逆变换、牛顿法 🧮

在本节课中，我们将学习如何对线性及非线性变换进行求逆，并深入探讨牛顿法这一强大的数值求解工具。我们将从线性变换的逆矩阵开始，逐步过渡到使用牛顿法求解非线性方程的根，最终实现复杂非线性变换的求逆。

线性变换与逆矩阵 🔄

上一节我们探讨了线性变换及其矩阵表示。本节中我们来看看如何“撤销”一个线性变换，即求其逆变换。

对于一个线性变换，如果存在另一个变换能将其效果完全抵消，使物体回到原始状态，那么这个变换就是可逆的，我们称后者为前者的逆变换。

在矩阵表示中，若矩阵 A 表示的变换是可逆的，则存在一个矩阵 A⁻¹，满足：
A⁻¹ * A = I 且 A * A⁻¹ = I
其中 I 是单位矩阵。

在Julia中，我们可以直接计算矩阵的逆：

A = [1 2; 3 4]
invA = inv(A)
# 验证 A * invA 是否近似等于单位矩阵
A * invA ≈ I

非线性变换的挑战与思路 🌀

然而，现实世界中的许多变换，如图像的透视变形、非线性扭曲等，都不是线性的。它们的函数形式复杂，通常无法简单地用一个矩阵来表示，因此也无法直接使用矩阵求逆的方法。

对于非线性函数 y = f(x)，求逆意味着：给定输出 y，我们要找到对应的输入 x，使得 f(x) = y。这等价于求解方程 f(x) - y = 0 的根。

牛顿法：一维情况 📉

为了求解非线性方程的根，我们引入牛顿法。这是一种迭代算法，通过线性逼近来快速收敛到方程的根。

牛顿法的核心思想是：从一个初始猜测值 x₀ 开始，利用函数在该点的切线（线性近似）来找到与x轴的交点，以此作为下一个、更接近真实根的估计值 x₁。

其迭代公式为：
xₙ₊₁ = xₙ - f(xₙ) / f'(xₙ)
其中 f'(xₙ) 是函数在 xₙ 处的导数。

以下是牛顿法在Julia中的简单实现，用于寻找函数 f(x) = x² - 2 的根（即 √2）：

function newton1d(f, fprime, x0; tol=1e-10, maxiter=100)
    x = x0
    for i in 1:maxiter
        fx = f(x)
        abs(fx) < tol && break # 如果f(x)足够接近0，则停止
        x = x - fx / fprime(x) # 牛顿迭代更新
    end
    return x
end

f(x) = x^2 - 2
fprime(x) = 2x # 导数
root = newton1d(f, fprime, 1.5) # 从1.5开始寻找根

牛顿法：扩展到多维与求逆 🧭

对于从二维到二维的非线性变换 y = T(x)，求逆问题同样可以转化为求根问题：给定 y，寻找 x 使得 G(x) = T(x) - y = 0。

此时，牛顿法的公式需要从标量形式推广到向量形式。迭代步骤变为：
xₙ₊₁ = xₙ - J⁻¹(xₙ) * G(xₙ)
其中 J(xₙ) 是变换 T 在点 xₙ 处的雅可比矩阵（即导数在多维的推广）。J⁻¹(xₙ) * G(xₙ) 在计算上等价于求解线性方程组 J(xₙ) * Δ = -G(xₙ) 得到 Δ，然后更新 xₙ₊₁ = xₙ + Δ。

在Julia中，我们可以利用自动微分来计算雅可比矩阵，并利用反斜杠运算符 \ 来高效求解线性方程组：

using ForwardDiff

function newton_inverse(T, y, x0; tol=1e-10, maxiter=100)
    x = x0
    for i in 1:maxiter
        G(x) = T(x) - y # 定义目标函数
        g_val = G(x)
        norm(g_val) < tol && break # 如果G(x)的范数足够小，则停止
        J = ForwardDiff.jacobian(G, x) # 自动计算雅可比矩阵
        delta = -J \ g_val # 求解 J * delta = -G(x)
        x = x + delta
    end
    return x
end

图像变换求逆的实践 🖼️

在图像处理中，我们经常需要根据变换后的像素位置，反推它在原图中的位置来获取颜色值（即逆向映射），以避免出现空洞或重叠。这正是非线性变换求逆的一个典型应用。

以下是实现图像非线性变换求逆的关键步骤概述：

1. 定义变换：明确将原图坐标 (x, y) 映射到新坐标 (x', y') 的函数 T。
2. 目标求逆：对于变换后图像上的每个像素点 p'，我们需要找到原图中的点 p，使得 T(p) = p'。
3. 应用牛顿法：以 p' 作为目标值 y，以一个合理的初始猜测（例如，假设变换不大，可用 p' 本身作为初始值）开始，利用上述的 newton_inverse 函数求解 p。
4. 采样颜色：根据求得的原图坐标 p（可能是非整数），通过插值方法（如双线性插值）从原图中获取颜色，并填充到变换后图像的 p' 位置。

总结 📚

本节课中我们一起学习了变换求逆的核心概念与方法。

• 我们回顾了线性变换的逆可以通过求逆矩阵得到。
• 我们认识到非线性变换求逆的复杂性，并将其转化为求根问题。
• 我们深入学习了牛顿法，从一维求根公式开始，理解了其通过线性逼近（切线） 快速迭代收敛的原理。
• 我们将牛顿法成功推广到多维情形，用于求解向量值函数的根，从而实现了对复杂非线性变换的数值求逆。
• 最后，我们看到了这一方法在图像扭曲与恢复中的实际应用，理解了逆向映射的重要性。

牛顿法是一个强大而基础的工具，其思想——用可解的线性问题去逼近和解决非线性问题——在科学计算的许多领域都有广泛应用。

L6：动态规划与接缝裁剪 (Seam Carving) 🧠✂️

在本节课中，我们将要学习动态规划的基本概念，并探索一个非常有趣的应用——接缝裁剪（Seam Carving）。这是一种智能的图像缩放技术，能够在不扭曲图像主要内容的前提下改变其尺寸。

概述

动态规划是一种解决优化问题的强大算法思想，其核心在于避免重复计算子问题。本节课我们将通过一个棋盘路径最小和的例子来理解动态规划，然后将其应用于接缝裁剪算法，学习如何智能地压缩图像。

动态规划简介

动态规划中的“规划”一词指的是解决优化问题，而非编写代码。它是一种算法设计范式，通过将复杂问题分解为重叠的子问题，并存储子问题的解来避免重复计算，从而高效地找到最优解。

路径最小和问题

为了理解动态规划，我们首先来看一个具体问题：在一个由随机数字（0-9）组成的 N x N 网格（如棋盘）中，寻找一条从顶部到底部的路径，使得路径上经过的数字之和最小。

路径的规则是：从顶部开始，每一步只能向正下方、左下方或右下方移动一格。这类似于国际象棋中“后”的移动方式，但只能向下。

对于一个 8x8 的网格，符合规则的路径总数是 11,814 条。一种简单（但低效）的方法是遍历所有路径并计算其和，然后找出最小值。

公式表示：
假设网格为矩阵 M，路径 P 是一系列坐标 (i, j) 的集合，其中 i 从 1 到 N，且 j 的相邻行间变化不超过 1。目标是找到路径 P*，使得：
sum(M[i, j] for (i, j) in P*) 最小。

重叠子问题与高效解法

上一节我们介绍了问题的定义和朴素解法。本节中我们来看看为什么动态规划更高效。

关键在于重叠子问题。假设我们想计算经过某个特定格子 (i, j) 的最小路径和。从该格子出发到底部的所有路径，会与从它下方三个相邻格子出发的路径大量重复。

动态规划思路：
我们不再从顶部开始向前寻找所有路径，而是从底部开始反向计算。

1. 对于最后一行，从每个格子出发的最小和就是该格子自身的值。
2. 对于上一行的每个格子，其最小路径和等于该格子的值加上其下方三个相邻格子中已计算出的最小路径和的最小值。
3. 依此类推，直到顶部。最终，顶部每个格子的值就代表了从该格子出发的最小路径和。全局最小值就是顶部这些值中的最小值。

这种方法将计算复杂度从指数级（检查所有路径）降低到了与网格像素数量成线性关系。

代码描述（伪代码）：

# 假设 energy 是一个 M x N 的矩阵，存储每个格子的“能量值”（在此问题中即格子本身的数字）
# 初始化一个同样大小的矩阵 dp 来存储最小路径和
dp = zeros(size(energy))
dp[end, :] = energy[end, :] # 最后一行直接复制

# 从倒数第二行开始向上迭代
for i in (M-1):-1:1
    for j in 1:N
        # 考虑下方相邻的列索引，注意边界处理
        left = max(1, j-1)
        right = min(N, j+1)
        # 找到下方相邻格子 dp 值的最小值
        min_below = minimum(dp[i+1, left:right])
        # 当前格子的 dp 值 = 自身能量 + 下方最小 dp 值
        dp[i, j] = energy[i, j] + min_below
    end
end
# 全局最小路径和就是 dp 第一行中的最小值
min_total_energy = minimum(dp[1, :])

接缝裁剪算法应用 🖼️

理解了动态规划后，我们来看一个精彩的应用——接缝裁剪。这是一种内容感知的图像缩放技术，目标是在改变图像宽高比时，尽可能保留重要的视觉内容。

算法核心思想

接缝裁剪的核心是定义图像中每个像素的“重要性”或“能量”。通常，边缘像素（颜色、亮度变化剧烈的区域）更为重要，因为它们通常对应物体的轮廓。

能量计算：
常用方法是计算图像的梯度幅值。梯度表示像素值变化最快的方向和速率，其幅值越大，说明该像素处于边缘的可能性越高，也就越重要。

1. 使用 Sobel 滤波器等算子计算图像在 x 和 y 方向上的梯度（近似偏导数）。
2. 每个像素的能量 E(i, j) 为其梯度幅值：E(i, j) = sqrt(Gx² + Gy²)。

寻找接缝：
一条“接缝”是从图像顶部到底部的一条连通路径，每行有且仅有一个像素，且相邻行的像素在列索引上相差不超过 1（即同样只能向下、左下或右下移动）。
我们的目标是找到一条能量总和最小的接缝。这正好转化成了我们之前解决的动态规划问题，其中网格是图像，每个格子的“值”就是其能量 E(i, j)。

移除接缝：
找到能量最小的接缝后，将其从图像中移除（将所有右侧的像素向左移动一列），图像宽度便减少了一像素。重复此过程，即可将图像缩放到任意宽度，同时尽量避开重要的边缘区域。

算法效果与局限性

接缝裁剪算法对于背景简单、主体突出的图像效果非常好（例如课程中演示的 MIT Stata 中心大楼和萨尔瓦多·达利的《记忆的永恒》）。它能巧妙地穿过天空、地面等平滑区域，避开建筑、钟表等轮廓清晰的物体。

然而，对于充满复杂纹理的图像（如文森特·梵高的《星夜》），由于处处都是“边缘”，算法可能无法找到明显的低能量路径，导致裁剪效果不理想。这时，可能需要设计更复杂的能量函数。

总结

本节课我们一起学习了动态规划的基本原理及其在接缝裁剪算法中的应用。

1. 动态规划通过解决重叠子问题来高效求解最优路径问题，避免了穷举法带来的指数级复杂度。
2. 接缝裁剪利用动态规划寻找图像中能量最低的连通路径（接缝）并移除它，从而实现智能的图像缩放。其关键在于使用梯度幅值来定义像素的能量，以保护重要的边缘信息。

这种将数学概念（梯度、优化）与算法思想（动态规划）相结合，解决实际图像处理问题的方法，正是计算思维的生动体现。

L7：结构 🧩

在本节课中，我们将学习如何识别和利用数据中的“结构”。结构意味着数据并非完全随机，而是存在某种模式或规律，我们可以利用这些规律来更高效地存储、处理和理解数据。我们将通过几个具体的例子来探索这个概念，包括稀疏矩阵、随机向量的统计特性以及主成分分析。

---

结构示例：从存储效率说起

上一节我们提到了利用结构可以优化算法。本节中，我们来看看如何利用结构来优化数据的存储。

独热向量

一个“独热向量”是一个向量，其中只有一个元素是1，其余所有元素都是0。例如：
[0, 1, 0, 0, 0, 0]

虽然这个向量有6个元素，但我们真的需要存储6个数字吗？实际上，我们只需要两个信息：向量的长度 n 和那个“1”所在的位置 k。在Julia中，我们可以创建一个自定义类型来利用这种结构。

代码示例：定义一个独热向量类型

struct OneHot <: AbstractVector{Int}
    n::Int  # 向量长度
    k::Int  # 非零元素位置
end

通过定义 size 和 getindex 方法，我们可以让这个 OneHot 类型表现得像一个普通向量，但内部只存储两个整数。

代码示例：让 OneHot 类型像向量一样工作

Base.size(v::OneHot) = (v.n,)
Base.getindex(v::OneHot, i::Int) = Int(i == v.k)

这样，当我们索引 OneHot(6, 2)[4] 时，会返回0；索引 OneHot(6, 2)[2] 时，会返回1。我们成功利用了结构，避免了存储所有零。

对角矩阵

对角矩阵是另一种具有明显结构的矩阵，只有对角线上的元素可能非零。对于一个 n x n 的对角矩阵，我们只需要存储 n 个数字，而不是 n^2 个。

Julia内置了 Diagonal 类型来处理这种结构。与存储所有元素的“稠密”矩阵相比，Diagonal 类型只存储对角线元素，节省了大量空间。

---

结构示例：数据中的模式

并非所有结构都是为了节省存储空间。有时，结构体现在数据本身的内在关系上。

随机向量的统计结构

一个由1到9的随机整数组成的百万维向量，看似完全随机，没有结构。但实际上，它拥有统计意义上的结构。

例如，我们可以计算它的均值和标准差。对于在1到9上均匀分布的随机变量，其理论均值是5，理论标准差是 sqrt((9^2 - 1)/12) ≈ 2.58。多次生成这样的随机向量，其样本均值和标准差都会稳定在这些理论值附近。

公式：均值与标准差

• 均值：mean(v) = sum(v) / length(v)
• 标准差：std(v) = sqrt( sum( (v .- mean(v)).^2 ) / (length(v) - 1) )

这种稳定性就是一种“结构”。我们有时甚至可以用均值和标准差这两个摘要统计量来代表整个数据集，而忽略具体的百万个数字。

乘法表（外积）

一个矩阵如果是由两个向量的外积构成的，即矩阵的每个元素 A[i, j] = v[i] * w[j]，那么这个矩阵就具有极强的结构。它被称为秩1矩阵。

对于一个 m x n 的秩1矩阵，我们只需要存储 m + n 个数字（向量 v 和 w），而不是 m * n 个。然而，仅凭观察矩阵的值，我们很难一眼看出它是否具有这种结构。

代码示例：检查并分解乘法表

function factor_multiplication_table(A)
    v = A[:, 1]  # 取第一列
    w = A[1, :]  # 取第一行
    scale = v[1]
    v_normalized = v ./ scale
    w_scaled = w .* scale
    # 检查重构矩阵是否近似等于原矩阵
    if isapprox(outer(v_normalized, w_scaled), A)
        return v_normalized, w_scaled
    else
        error("输入矩阵不是乘法表")
    end
end

---

揭示隐藏结构：奇异值分解与主成分分析

如果矩阵不是完美的秩1矩阵，而是“近似”的秩1矩阵（例如一个秩1矩阵加上一些随机噪声），或者是由多个秩1矩阵相加而成，我们如何发现其中的结构呢？这就需要用到强大的奇异值分解工具。

奇异值分解

奇异值分解可以将任意矩阵 A 分解为多个秩1矩阵（外积）的和：
A = U * Σ * V^T
其中，Σ 是一个对角矩阵，其对角线上的值（奇异值）从大到小排列。U 和 V 的列向量分别代表了数据在行和列方向上的主要模式。

最大的奇异值对应的 U 的第一列和 V 的第一列相乘，就给出了对原矩阵 A 的“最佳”秩1近似。前k个奇异值对应的部分相加，就给出了对 A 的最佳秩k近似。

主成分分析：数据的视角

主成分分析是SVD在统计学中的叫法，它为我们提供了一种理解高维数据的视角。

假设我们有一个数据矩阵，每一行是一个数据样本，每一列是一个特征。PCA的目标是找到数据中方差最大的方向（即主成分）。这等价于对数据中心化后的矩阵进行SVD。

可视化理解：

• 将一个秩1矩阵的两行数据分别作为x轴和y轴的值绘制成散点图，所有点会精确地落在一条直线上。
• 如果对这个矩阵添加一些噪声，数据点会散布在这条直线周围。
• PCA可以帮助我们找到这条“隐藏”的直线，即数据的主要变化方向。第一个主成分就指向这个方向。数据点沿主成分方向的散布（方差）很大，而垂直于主成分方向的散布（噪声）很小。

通过只保留前几个主成分，我们可以用更少的维度来近似表示原始数据，实现降维和去噪。这在图像压缩（如JPEG）、数据可视化以及许多机器学习任务中都有广泛应用。

---

总结

本节课中我们一起学习了“结构”在计算思维中的核心地位。

1. 存储结构：我们看到了如何利用独热向量、对角矩阵、稀疏矩阵的特定模式，来设计高效的数据存储方式，避免保存不必要的零值或冗余信息。
2. 数据模式：我们认识到即使是随机数据，也存在如统计摘要（均值、方差）这样的稳定结构。此外，乘法表（外积）是一种特殊的矩阵结构。
3. 结构发现：我们引入了强大的奇异值分解工具及其统计学版本——主成分分析。SVD/PCA能够自动发现数据中隐藏的层级结构（秩1分量），并允许我们根据重要性（奇异值大小）来近似表示数据，广泛应用于压缩、去噪和洞察数据本质。

理解并利用结构，是编写高效算法、构建智能模型和从数据中提取洞见的关键一步。

L8：主成分分析 📊

在本节课中，我们将学习如何将图像视为数据，并引入一种强大的数据分析工具——主成分分析。我们将从矩阵的角度理解数据，并利用线性代数中的变换工具来分析和提取数据中的关键结构。

上一节我们讨论了图像作为数据的表示方法，本节中我们来看看如何从统计和概率的视角分析这些数据。这是一个很好的过渡话题，我们将把矩阵视为数据，并尝试用我们已经见过的变换工具来分析它。

从图像到数据点

我们首先回顾一下之前的工作。我们希望通过某种方式理解数据。我们观察一个矩阵，并在作业中也会处理类似的问题。

我们之前讨论了外积。外积本质上是一个乘法表。例如，矩阵的第一列包含数字1到10，第二列则是第一列数字的两倍。这意味着矩阵的每一列都包含相同的信息，只是按不同倍数缩放。同样，每一行也是前一行的倍数。

这种矩阵被称为秩为一的矩阵。这意味着矩阵中本质上只包含一个向量的信息。虽然我们需要两个向量来存储它（一个列向量和一个缩放系数向量），但信息量是单一的。

在作业中，你将创建一个新的Julia类型来实现这种矩阵。该类型内部只存储两个向量，但在索引时表现得像一个完整的矩阵。这是抽象数组的一个很好的例子。

当我们可视化秩为一的矩阵时，会看到棋盘格图案，因为每一列看起来都相同，只是颜色强度根据缩放系数变化。当我们叠加另一个秩为一的矩阵（形成秩为二的矩阵）时，棋盘格图案开始消失。

在上一讲的最后，Alan介绍了奇异值分解这个来自线性代数和数值代数的工具。SVD可以将任何矩阵分解为多个秩为一的矩阵之和，从而帮助我们提取矩阵中的结构信息。

将图像行视为数据坐标

现在，我们开始从图像分析过渡到更通用的数据分析。我们将图像的行视为数据的坐标。

具体做法是：取图像的一个子块，将第一行作为X坐标，第二行作为Y坐标，然后将这些点绘制在普通的二维图表上。图表上的每个点对应原始图像矩阵的一列，该列包含两个数据点（X和Y坐标）。

我们绘制了两组点：

• 红色点：来自单个秩为一的矩阵（一个乘法表）。
• 蓝色点：来自同一个矩阵，但为每个元素添加了噪声（随机数）。

蓝色点云明显分散在一条直线周围。假设我们不知道这些数据来自一条带噪声的直线，我们的目标是重建这条直线。这就是统计学和机器学习中的核心问题：找到以某种标准衡量的“最佳”拟合直线。

一种常见的方法是最小二乘法，它最小化每个数据点到直线的垂直距离的平方和。然而，我们将采用一种略有不同的方法：最小化数据点到直线的垂直距离（即点到直线的最短距离）。这种方法能更好地泛化。

分析数据点云：均值与方差

为什么红色点会形成一条直线？因为对于原始图像中的每一列 (x, y)，y/x 是一个常数。所有点都位于一条通过原点的直线上。

对于蓝色的噪声点云，我们首先需要衡量它的大小。一个直观的想法是测量点云的“宽度”。但简单地取最大值和最小值之差容易受到异常值的影响。因此，我们采用统计方法。

首先，找到点云的中心，即计算均值。在Julia中，可以使用 Statistics 标准库中的 mean 函数。

计算出均值 (mean_x, mean_y) 后，我们将每个数据点减去均值，使点云中心移动到坐标原点 (0,0)。

现在，我们想测量中心化后数据在X方向上的宽度。由于数据有正有负，直接求平均会相互抵消得到0。我们需要一个衡量平均距离的量。

一种方法是计算平均绝对偏差，即求所有X坐标绝对值的平均值。这相当于将点云沿Y轴折叠。我们计算得到大约0.25。

然而，更常用的度量是标准差。其计算步骤如下：

1. 计算所有X坐标的平方。
2. 求这些平方值的平均值，得到方差。
3. 对方差取平方根，得到标准差 (σ)。

标准差衡量了数据点偏离均值的典型距离。我们通常在图表上绘制 ±2σ 的线，对于正态分布数据，大约95%的数据点会落在这个范围内。

捕捉数据方向：相关性

仅用X和Y方向的标准差（即一个“盒子”）来描述数据是不够的，因为它无法捕捉数据点在这个盒子内的排列方式。数据可能沿着从左下到右上的对角线分布，这意味着两个变量（例如两次考试成绩）是正相关的。

为了捕捉这种方向信息，我们需要旋转视角。我们的目标是找到数据分布最分散的方向（即方差最大的方向）和最集中的方向。

我们可以通过以下两种等价方式实现：

1. 旋转坐标轴：让新坐标轴的方向沿着数据点云最长的方向，然后计算数据点在新轴方向上的投影的方差。
2. 旋转数据：固定坐标轴，将数据点云进行旋转，然后计算每次旋转后数据在水平方向（X轴）的方差。

当我们系统性地改变旋转角度θ并计算每个方向上的方差时，会得到一个关于θ的函数图。该图呈正弦波形，存在最大值和最小值，分别对应数据分布最长和最短的方向。

通过计算（例如使用Julia的 argmax 函数），我们可以找到这些关键方向及其对应的方差值。这些方向被称为主成分，最长的方向是第一主成分。

主成分分析的几何意义

主成分分析的目标就是找到这些关键方向。为什么它很重要？以两次考试成绩为例，第一主成分（长轴方向）可能代表学生的“综合能力”得分，这个单一数字很大程度上能区分不同学生。而垂直于它的方向（短轴）代表两次考试的相对表现差异，可能包含的信息较少。

在更高维度（例如三维），一个秩为一的矩阵的所有点会落在一条通过原点的直线上，添加噪声后形成一个圆柱形的点云。一个秩为二的矩阵的点则分布在一个“厚平板”上。

奇异值分解：计算的利器

之前我们通过旋转和搜索来计算主成分，这种方法计算量较大。奇异值分解 提供了一种更优雅、更高效的计算方法。

SVD指出，任何矩阵 M 都可以分解为三个矩阵的乘积：
M = U * Σ * V^T
其中：

• U 和 V 是正交矩阵（代表旋转或反射）。
• Σ 是对角矩阵，其对角线上的元素称为奇异值（代表拉伸因子）。

从几何变换的角度看，这意味着任何线性变换都可以分解为：先旋转 (V^T)，再沿坐标轴拉伸 (Σ)，最后再旋转 (U)。对于二维数据，这相当于将一个圆变换成一个椭圆，SVD直接给出了这个椭圆的长短轴长度（奇异值）和方向（U 的列向量）。

我们通过一个交互示例演示了这一点：从一个单位圆盘内的随机点开始，应用矩阵 [2 1; 1 1] 定义的变换。SVD分解出的奇异值（如2.6和0.38）正好对应了变换后椭圆的长轴和短轴长度，而 U 矩阵的列向量则指示了这些轴的方向。

高维视角与总结

最后，我们将SVD应用于原始的 2×300 图像数据矩阵。这个矩阵可以理解为二维空间中的300个点，也可以理解为300维空间中的2个点。SVD中的 V 矩阵代表了一个在300维空间中的旋转。虽然无法直观可视化300维空间，但我们可以通过投影观察这种旋转在二维平面上的效果，就像观察夜空中旋转的星轨。

本节课中我们一起学习了：

1. 如何将图像数据转化为可分析的坐标点。
2. 使用均值、方差和标准差描述数据点云的中心和离散程度。
3. 认识到需要分析数据的方向性，引入了主成分分析的思想。
4. 通过奇异值分解这一强大工具，高效地计算数据的主成分和拉伸因子。
5. 理解了PCA在降维和提取数据关键特征中的核心作用。

从下节课开始，我们将更深入地探讨随机性、概率论以及如何将其应用于数据分析。

L9：抽样与随机变量 🎲

在本节课中，我们将开始更详细地探讨概率与随机性。我们将学习什么是随机变量，如何在Julia中生成随机数，以及如何通过可视化和统计方法来理解随机过程的行为。

---

随机数生成

上一节我们介绍了课程概述，本节中我们来看看如何在Julia中生成随机对象。Julia提供了一个非常通用的函数 rand 来生成随机性。

以下是 rand 函数的一些基本用法：

• rand(1:6)：从范围 1:6 中均匀地随机选取一个整数（模拟掷骰子）。
• rand([2, 3, 5, 7, 11])：从给定的数组中随机选取一个元素。
• rand("MIT")：从字符串中随机选取一个字符。
• rand('a':'z')：从字符范围 'a' 到 'z' 中随机选取一个字符。
• rand()：生成一个在 [0, 1) 区间内均匀分布的随机浮点数。
• rand(1:6, 10)：生成一个包含10个随机整数的数组，每个整数都来自范围 1:6。
• rand(distinguishable_colors(10), 5, 5)：生成一个5x5的矩阵，每个元素是从10种可区分颜色中随机选取的一种。

rand 函数非常灵活，它有很多不同的“方法”来处理不同类型的输入集合。

---

计数与概率

当我们多次重复一个随机实验（如抛硬币）时，我们常常需要统计每个结果出现的次数。我们可以使用 countmap 函数（来自 StatsBase 包）来轻松实现这一点。

例如，模拟抛10次硬币并统计正反面次数：

using StatsBase
tosses = rand(["heads", "tails"], 10)
toss_counts = countmap(tosses)

countmap 返回一个字典（Dict），其中键是可能的结果，值是该结果出现的次数。我们可以通过将次数除以总实验次数来估算每个结果的概率。

随着实验次数（如抛硬币次数）的增加，观测到的概率会趋近于其理论值（如0.5），这体现了“大数定律”。

---

非均匀抽样与伯努利试验

到目前为止，我们进行的都是均匀抽样，即每个结果被选中的概率相同。但有时我们需要模拟有偏的过程，例如一枚正面朝上概率为0.7的硬币。

以下是实现加权抽样的方法：

• 一种直观的方法是生成一个1到10的随机整数，如果数字≤7则返回“正面”，否则返回“反面”。
• 更通用的方法是利用 rand() 生成一个 [0,1) 的随机数，如果这个数小于给定的概率 p，则返回 true（代表成功），否则返回 false。这被称为以概率 p 进行的伯努利试验。

我们可以用简洁的Julia函数定义伯努利试验：

bernoulli(p) = rand() < p

运行这个函数多次，并统计 true 出现的比例，该比例会趋近于我们设定的概率 p。

---

随机变量与分布

一个随机变量是一个其值由随机现象结果决定的对象。我们关心的是它取每个可能值的概率。所有可能值及其对应概率的对应关系，称为该随机变量的概率分布。

为了直观理解分布，我们可以对实验数据进行可视化。

---

可视化：条形图与直方图

我们可以使用条形图来可视化分类数据（如骰子的点数）。

以下是绘制条形图的步骤：

1. 进行多次实验（如掷骰子）。
2. 使用 countmap 统计每个点数出现的次数。
3. 使用 bar 函数绘制条形图，其中x轴是类别（点数），y轴是对应的计数。

随着实验次数增加，每个条形的高度会趋近于期望的均匀分布。

当我们处理数值数据（如多个骰子的点数之和）时，我们使用直方图。直方图的条形宽度代表一个数值区间（“分箱”），条形面积（高度×宽度）与该区间内数据点的数量成正比。

使用 histogram 函数可以绘制直方图。通过设置 norm=true，可以归一化直方图，使得所有条形的总面积之和为1，这样y轴就表示概率密度。

---

中心极限定理

一个令人着迷的现象是：当我们对多个独立的随机变量（如多个骰子的点数）求和时，其总和的分布会呈现出特定的形状。

以下是观察这一现象的步骤：

1. 定义一个函数 roll_dice(n)，用于返回 n 个骰子点数的和。
2. 多次运行该实验，收集这些和值。
3. 绘制这些和值的直方图。

随着骰子数量 n 的增加，原始和值的分布会向右移动并展宽。为了看到收敛的形状，我们需要进行标准化：

• 中心化：从每个和值中减去其均值（data .- mean(data)）。
• 缩放：将中心化后的数据除以其标准差（sigma）。

经过 (data .- mean(data)) ./ std(data) 标准化处理后，无论 n 是多少，标准化后数据的分布都会收敛到同一个形状——著名的钟形曲线，即正态分布（或高斯分布）。

其概率密度函数公式为：

f(x) = (1 / √(2π)) * exp(-x² / 2)

这个现象就是中心极限定理的核心内容：大量独立随机变量之和的标准化形式，在极限条件下趋近于正态分布。这解释了为何正态分布在自然界和统计学中如此普遍。

---

本节课中我们一起学习了随机性的基本概念。我们探索了如何在Julia中生成随机数，如何进行计数和估算概率，以及如何实现非均匀抽样。我们定义了随机变量及其分布，并利用条形图和直方图进行可视化。最后，我们通过模拟多个骰子求和的实验，直观地验证了中心极限定理，观察到独立随机变量之和的标准化分布如何收敛到正态分布。这些工具和概念是理解更复杂随机过程和进行数据分析的基础。

L10：随机模拟建模 🎲

在本节课中，我们将学习如何使用概率来模拟随时间变化的系统。我们将建立一个简单的模型，用于模拟机械部件或灯泡的失效过程。通过这个过程，我们将探索如何定义新类型、扩展函数以及进行随机模拟。

---

概述 📋

我们将建立一个“灯泡失效”模型。在这个模型中，一组灯泡最初是完好的（绿色），在每个时间步，每个灯泡都有一定的概率失效（变为红色），并在下一个时间步永久失效（变为紫色）。我们将通过这个简单的随机过程，学习如何连接微观的个体行为与宏观的系统行为。

上一节我们介绍了随机模拟的基本概念，本节中我们来看看如何具体实现一个简单的随机过程模型。

---

模型描述与模拟 🔧

我们有一个灯泡阵列，每个灯泡的状态随时间变化。模型的核心规则如下：

• 初始时，所有灯泡都是完好的（状态为0）。
• 在每个时间步 t，每个仍完好的灯泡以概率 p 失效。
• 如果一个灯泡在时间 t 失效，我们将其状态标记为 t（表示失效时间）。
• 失效后的灯泡在后续时间步状态保持不变（永久失效）。

以下是模拟这个过程的Julia代码核心逻辑：

function simulate_failure(N, p, max_time)
    # 初始化状态矩阵，0表示完好
    state = zeros(Int, N, N)
    for t in 1:max_time
        for i in 1:N, j in 1:N
            # 如果灯泡完好，则进行随机测试
            if state[i, j] == 0 && rand() < p
                state[i, j] = t  # 记录失效时间
            end
        end
    end
    return state
end

---

可视化绘图 🎨

为了直观展示模拟结果，我们需要绘制灯泡阵列。这涉及到在特定坐标上绘制形状（如圆形和方形）。

以下是绘制矩形和圆形的辅助函数：

function rectangle(x, y, w, h)
    xs = [x, x+w, x+w, x]
    ys = [y, y, y+h, y+h]
    return (xs, ys)
end

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_69.png)

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_71.png)

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_73.png)

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_75.png)

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_77.png)

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_79.png)

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_81.png)

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_83.png)

![](https://github.com/OpenDocCN/cs-notes-zh/raw/master/docs/mit-18s191-intro/img/b8e1d125d9ff89215ae2d0f3136f14a4_85.png)

function circle(x, y, r, n=30)
    θ = range(0, 2π, length=n)
    xs = x .+ r .* cos.(θ)
    ys = y .+ r .* sin.(θ)
    return (xs, ys)
end

在绘图时，我们可以使用 fill=true 参数来填充形状，并使用 annotate! 函数添加文本标签。字符串插值（$）可以方便地将变量值嵌入标签字符串中。

---

数据分析与宏观视角 📊

运行模拟后，我们可以统计每个时间步处于不同状态（绿色/红色/紫色）的灯泡数量。

以下是可能的数据分析步骤：

1. 统计数量：遍历状态矩阵，根据存储的值（0表示完好，t 表示在时间 t 失效）对灯泡进行分类计数。
2. 绘制曲线：将完好灯泡数量、新失效灯泡数量随时间变化的曲线绘制出来。
3. 观察规律：我们会发现，完好灯泡数量的曲线相对平滑，而每个时间步新失效数量的曲线噪声较大。这是因为累计数量（完好总数）对随机波动有平滑作用。

单个模拟的结果具有随机性。为了得到更稳定的规律，我们可以进行多次独立模拟并计算平均值。

---

随机变量与类型定义 🎲

上述模型中的每个灯泡，其“是否失效”的行为可以用一个伯努利随机变量来建模。伯努利随机变量模拟了一次加权硬币抛掷，它以概率 p 取值为1（成功/失效），以概率 1-p 取值为0（失败/完好）。

在Julia中，我们可以创建一个新类型来封装这个随机变量及其参数。

struct Bernoulli
    p::Float64  # 成功概率
end

接着，我们可以为这个新类型扩展已有的通用函数。例如，定义如何从 Bernoulli 对象中随机采样：

import Base: rand
function rand(b::Bernoulli)
    return Int(rand() < b.p)
end

我们还可以为定义在 Statistics 包中的 mean 函数添加方法，来计算伯努利随机变量的理论均值：

using Statistics
function mean(b::Bernoulli)
    return b.p
end

这样，我们就创建了一个计算实体，它既包含了数据（参数 p），也包含了对其进行的操作（采样、计算均值）。

---

从微观到宏观：确定性方程 🔬

通过对大量独立灯泡的微观行为进行平均，我们可以推导出描述系统宏观行为的确定性方程。

设 I_t 为时间 t 时完好灯泡的数量。

• 在时间 t，有 I_t 个灯泡完好。
• 每个灯泡失效的概率为 p。
• 因此，在时间 t 失效的灯泡数量的期望值是 p * I_t。

由此，我们可以得到 I_t 随时间变化的递推关系：
I_{t+1} = I_t - p * I_t = (1 - p) * I_t

这是一个确定性方程。假设初始数量为 I_0，我们可以直接解出：
I_t = (1 - p)^t * I_0

这表明完好灯泡的数量随时间呈指数衰减。当我们将这个确定性方程的解与多次随机模拟的平均结果进行比较时，会发现它们完美吻合。这揭示了在大量个体情况下，随机系统的平均行为可以由一个简单的确定性模型来描述。

---

总结 🏁

本节课中我们一起学习了：

1. 构建随机过程模型：我们建立了灯泡失效的简单随机模型，其中每个个体在每个时间步以固定概率改变状态。
2. 实现模拟与可视化：我们编写了模拟代码，并学习了如何绘制自定义形状来可视化模拟结果。
3. 定义抽象类型：我们将伯努利随机变量定义为一种新的Julia类型，并为其扩展了 rand 和 mean 函数，体现了将数据与操作封装在一起的思想。
4. 连接微观与宏观：我们通过分析推导出，大量独立个体随机行为的平均效应，可以用一个描述总数量的确定性指数衰减方程来刻画。这是理解复杂系统的一种强大方法。

这个简单的模型是许多领域（如流行病学、金融、化学）中更复杂模型的基础。通过掌握这些基本概念，你便拥有了分析随机动态系统的起点。

L11：作为类型的随机变量 🎲

在本节课中，我们将学习如何利用Julia的类型系统来优雅地表示和操作随机变量。虽然内容以概率统计为背景，但其核心是关于如何通过抽象类型和子类型来组织代码，这是一种强大的计算机科学思想。

概述

我们将从随机变量的基本概念出发，逐步构建一个类型系统来表示高斯分布、伯努利分布等。你将看到，通过定义抽象类型和具体方法，我们可以轻松实现随机变量的命名、参数化、采样、求均值/方差以及相加等操作。这种设计模式不仅限于随机变量，可以应用于许多需要抽象和层次化组织的软件系统。

随机变量与概率分布

随机变量是概率论中的核心概念。对于一个离散随机变量，我们为其每个可能的结果分配一个概率。这组概率的集合被称为概率分布。

以高斯分布（正态分布）为例，它由两个参数决定：均值（μ）和标准差（σ）。均值决定了分布的中心位置，而标准差决定了分布的宽度。当μ=0，σ=1时，我们得到标准正态分布。

在Julia中，我们可以用以下方式生成标准正态分布的样本：

data = randn(10^5)

要生成均值为mu、标准差为sigma的样本，可以进行变换：

data = mu .+ sigma .* randn(10^5)

高斯分布的概率密度函数（PDF）公式为：

pdf(x) = (1 / (σ * √(2π))) * exp(-(x - μ)^2 / (2σ^2))

高斯分布的可加性

高斯分布有一个美妙的性质：两个独立高斯随机变量之和仍然是高斯随机变量。新分布的均值是原均值之和，方差是原方差之和。

如果 X ~ N(μ₁, σ₁²)， Y ~ N(μ₂, σ₂²)，且X与Y独立，那么：

X + Y ~ N(μ₁ + μ₂, σ₁² + σ₂²)

这个性质使得处理高斯分布的和变得非常简单。

设计随机变量的类型系统

上一节我们介绍了随机变量的基本操作，本节中我们来看看如何用类型系统来抽象地表示它们。关键在于“抽象”与“具体”的层次设计。

首先，我们定义一个抽象类型RandomVariable，作为所有随机变量的根。

abstract type RandomVariable end

然后，我们可以定义其子类型，例如区分连续型和离散型随机变量：

abstract type DiscreteRandomVariable <: RandomVariable end
abstract type ContinuousRandomVariable <: RandomVariable end

符号 <: 表示“是……的子类型”。

实现具体分布：高斯分布

现在，我们来实现第一个具体分布——高斯分布。它属于连续型随机变量。

定义一个高斯分布类型需要存储其均值和方差：

struct Gaussian <: ContinuousRandomVariable
    μ
    σ²
end

我们可以为它定义一个便捷的构造函数，创建标准正态分布：

Gaussian() = Gaussian(0, 1)

定义了类型，我们就可以为其实现各种方法。例如，计算理论均值和方法：

mean(d::Gaussian) = d.μ
var(d::Gaussian) = d.σ²

对于任何随机变量，标准差都是方差的平方根。我们可以为顶层的RandomVariable定义这个通用方法：

std(d::RandomVariable) = sqrt(var(d))

这是一个“泛化”的例子：我们为抽象类型定义通用行为，所有子类型只要实现了var方法，就能自动获得std方法。

随机变量的加法

我们希望能够将随机变量相加。对于高斯分布，我们可以利用其可加性，直接返回一个新的高斯分布：

import Base: +
+(X::Gaussian, Y::Gaussian) = Gaussian(mean(X) + mean(Y), var(X) + var(Y))

但是，并非所有随机变量相加后都能得到同类型的简单分布。例如，两个伯努利分布相加的结果不再是伯努利分布。

为了解决通用性问题，我们创建一个新的类型来表示“两个随机变量的和”：

struct SumOfTwoRandomVariables <: RandomVariable
    X
    Y
end

然后，我们可以为+运算符定义一个更通用的方法，当不知道具体规则时，就返回这个“和”类型：

+(X::RandomVariable, Y::RandomVariable) = SumOfTwoRandomVariables(X, Y)

对于“和”类型，我们也可以定义其均值等方法。根据期望的线性性质，无论X和Y是什么，和的均值等于均值之和：

mean(S::SumOfTwoRandomVariables) = mean(S.X) + mean(S.Y)

如果X和Y独立，方差也可加：

var(S::SumOfTwoRandomVariables) = var(S.X) + var(S.Y)

最后，为了能从“和”类型中采样，我们定义其rand方法：分别从两个组成部分采样然后相加。

rand(S::SumOfTwoRandomVariables) = rand(S.X) + rand(S.Y)

实现更多分布：伯努利分布

现在，我们添加另一个常见分布——伯努利分布。它是一个离散随机变量，只有一个参数p（成功概率）。

以下是伯努利分布的定义：

struct Bernoulli <: DiscreteRandomVariable
    p
end

其理论均值是p，方差是p*(1-p)：

mean(B::Bernoulli) = B.p
var(B::Bernoulli) = B.p * (1 - B.p)

从伯努利分布采样的方法是：生成一个[0,1)区间的均匀随机数，若小于p则返回1，否则返回0。

rand(B::Bernoulli) = rand() < B.p ? 1 : 0

类型系统的威力

有了上述基础，类型系统的威力开始显现。我们可以轻松创建复杂的随机变量表达式。

例如，创建一个高斯分布与一个伯努利分布的和：

mixture = Gaussian(4, 0.3) + Bernoulli(0.7)

我们可以计算这个混合分布的理论均值：

mean(mixture) # 等于 mean(Gaussian(4,0.3)) + mean(Bernoulli(0.7))

我们也可以从这个混合分布中采样并绘制直方图，观察其形状。由于我们定义了SumOfTwoRandomVariables的rand方法，这一切都能自动工作。

更强大的是，由于+运算符是通用的，我们可以轻松地添加多个随机变量：

sum_of_three = Bernoulli(0.25) + Bernoulli(0.25) + Gaussian(0, 1)

这实际上构建了一个表达式树（SumOfTwoRandomVariables内部嵌套另一个SumOfTwoRandomVariables），但所有关于均值、方差和采样的方法都能递归地正确工作。

我们甚至可以使用Julia内置的sum函数对一组随机变量求和：

sum([Bernoulli(0.25) for i in 1:30])

这等价于连续使用+运算符，得益于我们为+定义的通用方法。

扩展：卡方分布

为了展示该框架的扩展性，我们再添加一个分布：卡方分布。

卡方分布与高斯分布有关。自由度为1的卡方分布（记作χ²(1)）就是一个标准正态随机变量的平方。

struct Chisq1 <: ContinuousRandomVariable end

rand(::Chisq1) = randn()^2

自由度为k的卡方分布（χ²(k)）是k个独立的χ²(1)分布之和。利用我们已经建立的类型系统，我们可以轻松表示它：

# χ²(4) 分布
chisq4 = sum([Chisq1() for i in 1:4])

我们可以直接对chisq4进行采样和绘制直方图，而无需手动实现χ²(4)的复杂公式。这体现了通过组合简单部件来构建复杂系统的力量。

符号计算的可能性

我们构建的类型系统甚至能与符号计算结合。假设我们定义一些符号变量作为分布参数：

using Symbolics
@variables μ1, σ1², μ2, σ2²

然后，我们可以创建包含符号参数的高斯分布，并将它们相加：

G1 = Gaussian(μ1, σ1²)
G2 = Gaussian(μ2, σ2²)
G_sum = G1 + G2

此时，mean(G_sum)会正确地返回符号表达式μ1 + μ2，var(G_sum)会返回σ1² + σ2²。这表明，我们为数值计算设计的抽象，在符号计算中同样有效，这是泛型编程强大之处的体现。

总结

本节课中我们一起学习了如何将随机变量视为类型，并利用Julia的类型系统构建一个灵活而强大的概率建模框架。

我们首先定义了RandomVariable抽象类型及其子类型。然后为具体的Gaussian和Bernoulli分布实现了mean、var、rand等核心方法。通过定义SumOfTwoRandomVariables类型和通用的+运算符，我们实现了随机变量的加法，并自动处理了理论性质（如均值可加）和采样。

这种设计模式的核心优势在于：

1. 组织性：代码结构清晰，贴近数学概念。
2. 可扩展性：添加新分布（如Chisq1）非常容易，并能立即与现有系统（如加法）集成。
3. 泛化能力：在抽象类型上定义的方法（如std）自动适用于所有子类型。
4. 组合性：通过组合简单的分布可以构建复杂的分布，无需为每种复杂情况重写代码。

这堂课虽然以随机变量为例，但其展示的通过抽象类型和多重分派来构建可扩展、可组合的领域特定系统的思想，是Julia语言的核心优势，也能广泛应用于其他领域的软件设计。

L12：随机游走 1

在本节课中，我们将要学习随机游走的概念、应用及其在Julia中的实现方法。随机游走是模拟随机过程的基础模型，广泛应用于物理、金融和生物等领域。

什么是随机游走？

随机游走描述了一个粒子在空间中随机移动的过程。在每一步，粒子都会随机选择一个方向进行移动。这个过程可以用来模拟许多看似随机的现象。

上一节我们介绍了随机游走的基本概念，本节中我们来看看如何用Julia实现它。

随机游走的可视化

为了直观理解随机游走，我们可以观察一个在二维网格上移动的粒子。每一步，它都随机选择上、下、左、右四个方向之一移动。

随着步数的增加，粒子的轨迹会变得复杂，并且看起来越来越连续，仿佛在连续空间中运动。这个过程与气体分子的布朗运动非常相似。

随机游走的实际应用

随机游走模型在许多领域都有应用。例如，在金融中，股票价格的波动有时可以被建模为一种带有漂移的随机游走。在物理学中，它描述了微观粒子在流体中的扩散过程。在生物学中，它可以模拟中性基因在种群中的传播。

在Julia中实现随机游走

现在，我们来看看如何在Julia中实现一个简单的一维随机游走。最简单的模型是：粒子从0开始，每一步以1/2的概率向左（-1）或向右（+1）移动。

以下是生成单步随机移动的几种方法，我们将比较它们的性能。

# 方法1：使用元组
step1() = rand((-1, 1))

# 方法2：使用布尔运算
step2() = 2 * rand(Bool) - 1

# 方法3：使用比较
step3() = rand() < 0.5 ? -1 : 1

# 方法4：使用符号函数
step4() = sign(randn())

# 方法5：使用数组（不推荐，性能较差）
step5() = rand([-1, 1])

使用BenchmarkTools包可以方便地比较这些函数的性能。通常，使用元组或布尔运算的方法性能较好，而每次创建新数组的方法性能较差。

生成随机游走轨迹

有了单步移动函数，我们就可以生成整个随机游走的轨迹了。以下是一个生成一维随机游走轨迹的函数。

function walk_1d(n)
    x = 0
    xs = [x]
    for i in 1:n
        x += step1() # 这里可以使用上面定义的任何一种step函数
        push!(xs, x)
    end
    return xs
end

运行这个函数，并绘制多条轨迹，我们可以看到随机游走的扩散行为：轨迹会形成一个以原点为中心、随时间逐渐展宽的“云”。

迈向通用编程

上面的walk_1d函数是专门为一维情况编写的。如果我们想模拟二维或更高维的随机游走，就需要重写函数。为了避免代码重复，我们可以编写一个通用的随机游走函数。

通用函数的核心思想是将初始化和单步移动作为参数传入。

function generic_walk(initialize, step, n)
    x = initialize()
    trajectory = [x]
    for i in 1:n
        x += step()
        push!(trajectory, x)
    end
    return trajectory
end

这样，我们只需要定义不同维度的initialize和step函数，就可以用同一个generic_walk函数来生成轨迹。

使用类型构建随机游走器

为了进一步组织代码，我们可以定义表示随机游走器的类型。这让我们能将数据（位置）和行为（移动）封装在一起。

首先，定义一个抽象的Walker类型。

abstract type Walker end

然后，定义一维游走器类型及其相关方法。

struct Walker1D <: Walker
    pos::Int64
end

# 获取位置
position(w::Walker1D) = w.pos

# 定义移动步长（注意：这里返回的是位移量，不是新状态）
step(::Walker1D) = rand((-1, 1))

# 更新游走器状态（函数式风格，创建新对象）
function update(w::W, s) where {W <: Walker}
    new_pos = position(w) + s
    return W(new_pos)
end

使用where语法和类型参数W，update函数可以适用于任何Walker的子类型，并返回对应类型的新对象。

可变与不可变类型

在上面的例子中，Walker1D是不可变结构体。更新时，我们创建了一个新的游走器对象。另一种方法是使用可变结构体。

mutable struct MutableWalker1D <: Walker
    pos::Int64
end

# 原地更新函数（约定俗成，函数名以`!`结尾表示会修改参数）
function update!(w::MutableWalker1D, s)
    w.pos += s
    return nothing
end

使用可变对象在概念上更简单，但函数式风格（创建新对象）的代码通常更易于推理和调试，因为它避免了“副作用”。

通用的轨迹生成函数

最后，我们可以编写一个完全通用的轨迹生成函数，它接受任何Walker子类型的实例。

function trajectory(w::W, n) where {W <: Walker}
    pos_vec = [position(w)]
    current_walker = w
    for i in 1:n
        s = step(current_walker)        # 获取位移
        current_walker = update(current_walker, s) # 更新状态，得到新walker
        push!(pos_vec, position(current_walker))
    end
    return pos_vec
end

通过为不同的Walker子类型定义特定的position、step和update方法，这个通用的trajectory函数就能处理各种随机游走器。

总结

本节课中我们一起学习了随机游走的基本概念及其广泛的应用。我们探讨了在Julia中实现随机游走的多种方法，从简单的函数到通用的、基于类型的编程模式。我们比较了不同实现方式的性能，并介绍了通过定义类型和方法来实现通用、可扩展代码的技术。在下一节课中，我们将继续探索随机游走更深入的性质和应用。

L13：随机游动 2

在本节课中，我们将继续学习随机游动。我们将了解如何自定义Julia中对象的显示方式，探索累积求和函数 cumsum，学习如何创建和连接向量，并介绍两种三维绘图方法。

概述：自定义显示与结构化矩阵

上一节我们介绍了随机游动的基本概念。本节中，我们来看看如何通过自定义对象的显示方式来更好地理解数据结构，特别是结构化矩阵。

Julia允许我们自定义对象的显示方式。例如，我们可以定义一个类型，并指定其打印格式。

struct MyType
    value::Int
end

import Base.show
function Base.show(io::IO, obj::MyType)
    print(io, "这是MyType类型，值为 $(obj.value)")
end

y = MyType(10)  # 输出：这是MyType类型，值为 10

Julia为一些具有特定结构的矩阵（如下三角矩阵）预定义了类型。当显示这些矩阵时，零元素会以点（.）表示，使结构更清晰。

using LinearAlgebra
# 假设P是一个矩阵
L = LowerTriangular(P)  # 以更清晰的方式显示下三角部分

帕斯卡三角形与随机游动

帕斯卡三角形中的数字是二项式系数，可以通过 binomial(n, k) 函数计算。有趣的是，当我们只关注其中的奇数时，会呈现出谢尔宾斯基三角形的图案。

帕斯卡三角形的构建规则与随机游动有密切联系。每一行的新数字由上一行相邻两个数字相加得到，这类似于一个卷积操作。

随机游动中，物体在每一步以相等概率向左或向右移动。其位置随时间的变化是这些独立随机步长的累积和。

累积求和：cumsum 函数

在随机游动中，我们需要计算一系列随机步长的累积和，以得到运动轨迹。Julia中的 cumsum 函数正是用于此目的。

以下是 cumsum 函数的一个示例：

steps = [-1, 1, 1, -1, 1]
trajectory = cumsum(steps)  # 结果为：[-1, 0, 1, 0, 1]

为了得到从时间0开始的完整轨迹，我们可以在步长向量前连接一个0：

full_trajectory = [0; cumsum(steps)]

在Julia中，我们可以自己实现一个高效的 cumsum 函数，其性能与内置函数相当：

function my_cumsum(v)
    new_v = similar(v)
    new_v[1] = v[1]
    for i in 2:length(v)
        new_v[i] = new_v[i-1] + v[i]
    end
    return new_v
end

向量与轨迹的可视化

我们可以生成多条随机游动轨迹并进行可视化。每条轨迹本身就是一个随机对象，即一个随机函数。

当我们固定一个时间点（例如 t=25）并观察大量轨迹在该时刻的位置时，这些位置会形成一个概率分布。根据中心极限定理，这个分布会接近正态（高斯）分布。

概率演化的主方程

我们可以用数学公式描述随机游动位置概率分布随时间的演化。设 pᵗ 为一个向量，其中元素 pᵗᵢ 表示在时间 t 位于位置 i 的概率。

从时间 t 到 t+1 的演化遵循以下主方程：

pᵗ⁺¹ᵢ = (1/2) * pᵗᵢ₋₁ + (1/2) * pᵗᵢ₊₁

这个方程说明，在时间 t+1 到达位置 i 的概率，等于在时间 t 位于 i-1 并向右跳的概率，加上在时间 t 位于 i+1 并向左跳的概率（每次跳跃概率均为 1/2）。

这个演化过程本质上也是一个卷积操作。我们可以在Julia中实现这个演化过程，并观察概率分布如何随时间扩散。

function evolve(p)
    new_p = similar(p)
    n = length(p)
    # 处理边界（例如吸收边界）
    new_p[1] = 0.5 * p[2]
    new_p[end] = 0.5 * p[end-1]
    for i in 2:n-1
        new_p[i] = 0.5 * (p[i-1] + p[i+1])
    end
    return new_p
end

通过迭代这个方程，并从初始分布（例如，所有概率集中在原点）开始，我们可以模拟出概率波包的扩散过程，并将其在三维中可视化（时间 vs. 位置 vs. 概率）。

总结

本节课中我们一起学习了随机游动的更多内容。我们了解了如何利用Julia的类型系统自定义对象显示，探索了帕斯卡三角形与随机游动之间的深刻联系，并实践了使用 cumsum 函数计算轨迹。最重要的是，我们引入了描述随机游动概率分布演化的主方程，并通过模拟和可视化观察了概率的扩散过程。这些概念将随机游动从一个简单的路径模拟，提升到了对其整体统计性质的理解。

L14：离散与连续

在本节课中，我们将探讨离散数学与连续数学之间的相互作用。我们将看到，尽管它们常被视为不同的领域，但它们的结合往往能产生更强大的洞察力。我们将通过具体的例子，如计算圆周率面积和模拟随机游走，来理解离散与连续如何相互转化和补充。

离散与连续的定义

上一节我们介绍了课程的主题，本节中我们来看看如何定义“离散”与“连续”。在不涉及高深数学的情况下，我们通常通过例子来理解它们。

以下是离散数学对象的例子：

• 有限集合，例如序列 1 到 100。
• 无限离散集合，例如整数集。

相比之下，连续数学涉及极限、整个区间和曲面。例如，实数线包含了像 π 和 √3 这样的数。

离散与连续的融合

我们经常听到有人偏爱离散数学或连续数学。然而，在当今时代，这两个领域的界限正变得模糊。例如，机器学习将连续优化和梯度等概念引入了计算机科学。同样，数据科学和统计学也广泛应用连续数学思想。

连续数学有时比离散数学更简单。对于一个非常大的离散对象，将其转化为连续对象通常能保留核心特征并简化不必要的细节。

更重要的是，离散与连续的结合比单独使用任何一种都更有用。能够同时以两种方式思考会带来巨大的优势。在气候变化和流行病建模等现实应用中，连续数学也至关重要。

索引与函数求值

现在，让我们通过一个简单的例子来探讨离散与连续如何协同工作：索引和函数求值。

考虑一个向量 v 和取出其第 i 个元素的操作 v[i]。我们通常不将其视为函数，而只是“取出”一个元素。另一方面，对于函数 f(x)，如 sin(x)，我们将其视为一个“函数机器”，需要一个求值过程。

然而，从概念上讲，向量就是一个离散函数，其定义域是索引 1 到 n，函数值就是 v[i]。输入索引 i，输出元素 v[i]。

让我们用代码来具体说明。假设我们有一个包含数字 2 到 20 的向量。

v = [2:2:20;]  # 创建向量 [2, 4, 6, ..., 20]
v[7]  # 索引：取出第7个元素，结果是14

这里的 v[7] 感觉上是从内存中“取出”第7个元素。

但是，如果我们直接对一个“范围”对象进行索引：

r = 2:2:20    # 这是一个范围对象，只存储起始值、步长和结束值
r[7]          # 函数求值：计算 2 * 7，结果是14

这里的 r[7] 实际上是在计算函数 f(i) = 2i 在 i=7 时的值。这是一个明确的函数求值过程，尽管其定义域仅限于 1 到 10 的整数。相比之下，一个普通的函数 f(x) = 2x 可以接受任何实数输入，如 f(π)。

无论如何，向量在本质上就是一个函数。

圆的面积：离散近似与连续极限

上一节我们看到了索引如何隐含函数关系，本节中我们来看看如何用离散方法逼近连续对象的面积。

我们通过内接正多边形来逼近单位圆的面积（π）。随着多边形边数的增加，其面积越来越接近 π。

以下是不同边数多边形面积的近似值：

• 正方形（4边）：面积 ≈ 2.0
• 正八边形（8边）：面积 ≈ 2.828
• 正十六边形（16边）：面积 ≈ 3.061
• 正一百边形（100边）：面积 ≈ 3.139
• 正一千零二十四边形（1024边）：面积 ≈ 3.141

我们通常认为需要边数趋于“无穷”才能得到精确的 π，而“无穷”感觉非常遥远。然而，通过一种称为“外推”的技巧，我们可以从有限的数据中提取出更多关于极限的信息。

如果我们有一系列边数为 2^k（如 4, 8, 16, ..., 1024）的多边形面积值，我们可以对它们进行特定的线性组合（卷积）。例如，用 4/3 乘以当前面积减去 1/3 乘以前一个面积。

这样操作后，新的数列会比原始面积值更快地逼近 π，有效数字位数更多。我们可以重复这个过程，使用不同的系数（如 16/15, 64/63），每次都能揭示出更多隐藏的 π 的位数。

这个现象背后的原理是面积公式的泰勒展开。面积 A(s) 关于边数 s 的展开式包含 π 以及 s^(-2), s^(-4) 等项。当我们加倍边数时，s^(-2) 项变为原来的 1/4。精心选择的线性组合可以消去低阶误差项（如 s^(-2)），从而加速收敛。

这个例子表明，结合对连续极限（泰勒展开）的理解和离散计算，我们可以更高效地从离散数据中提取连续极限的信息。

面积的本质：多种离散途径通向同一连续极限

我们也可以用其他离散方法估算圆的面积，比如在圆上覆盖网格并计算内部的小方格数量。

随着方格越来越小，方格面积之和也会趋近于 π。这里有一个更深层的哲学或数学观点：如果我们暂时忘掉“面积”的直观概念，如何确信用不同离散方法（多边形、方格、甚至其他形状）得到的极限值是相同的？

逻辑上并不显然。数学上，正是由于许多不同的、合理的离散逼近过程都收敛到同一个数值，我们才确信这个共同的极限值是一个良定义的数学对象，并值得被称为“面积”。这种“多种途径通向同一终点”是定义许多连续数学概念（如积分、长度）的基础。

随机游走：从离散到连续的布朗运动

现在，让我们将“多路径收敛定义对象”的理念应用到随机过程中。

我们回顾一下帕斯卡三角形，它可以被归一化并绘制成一条曲线，表示在特定时间步，一个简单随机游走（掷硬币决定左右）所处位置的概率分布。

我们也可以模拟一个连续版本的随机游走。不是在每个时间步向左或向右移动固定距离，而是在每个微小的时间间隔内，根据一个均值为 0、方差与时间间隔成正比的正态分布来移动。

以下是模拟步骤：

1. 将总时间（例如 1）分为 N 小段。
2. 在每一小段，从正态分布 N(0, 1/N) 中抽取一个随机数作为位移。
3. 将这些位移累加，得到一条在时间上连续的路径。

当我们不断增加 N（细化时间间隔）时，这些离散生成的路径会趋近于一个连续的随机过程，称为布朗运动（或维纳过程）。

关键点在于：布朗运动的统计性质（如任意时刻的分布、任意两时刻的关联）是唯一的。无论底层离散步骤是基于正态分布还是伯努利分布（掷硬币），只要方差设置正确，最终的连续极限都是同一个布朗运动。

因此，布朗运动作为一个连续的数学对象“存在”，因为所有合理的离散逼近都导向它。在极限下，这个连续描述（用微分方程描述）通常比处理庞大的离散概率分布（如帕斯卡三角形）更简洁。

帕斯卡三角形的连续极限正是著名的钟形曲线（正态分布）。描述其随时间演化的离散规则（将相邻值相加并除以2）在连续极限下，就变成了描述扩散现象的偏微分方程（热方程）。

总结

本节课中我们一起学习了离散数学与连续数学之间的深刻联系。我们看到：

1. 离散与连续并非对立，而是相辅相成。
2. 连续极限通常能简化对大型离散系统的描述。
3. 通过结合离散计算和对连续极限的理解（如外推法），我们可以更高效地获取信息。
4. 许多连续的数学概念（如面积、布朗运动）之所以被良好定义，正是因为多种离散逼近方法都收敛到同一个结果。离散与连续是同一枚硬币的两面，掌握两者能提供更强大的问题解决视角。

L15：线性模型与模拟 📊

概述

在本节课中，我们将学习线性模型和统计模拟。我们将从处理带标签的数据开始，然后深入探讨如何为数据拟合一条直线，并理解统计输出表格中各项数字的含义。通过使用Julia进行快速模拟，我们将直观地理解线性回归背后的核心概念，例如系数估计、标准误差和假设检验。

---

数据组织：数据框 📋

上一节我们介绍了课程主题，本节中我们来看看如何组织数据。Julia借鉴了R语言，提供了数据框来处理带标签的数据。数据框的每一列可以有不同的数据类型，并且拥有列标签，这比普通的矩阵更便于处理现实世界的数据。

以下是如何创建一个数据框：

using DataFrames
# 假设 x 是华氏度数据，y 是摄氏度数据
data = DataFrame(degrees_fahrenheit = x, degrees_centigrade = y)

你可以轻松地在数据框和矩阵之间转换：

# 从矩阵创建数据框
df = DataFrame([x y], [:fahrenheit, :celsius])

# 将数据框转换回矩阵
matrix_data = Matrix(df)

数据框可以方便地读写到CSV文件：

using CSV
# 写入文件
CSV.write("test.csv", data)
# 从文件读取
data_from_file = CSV.read("test.csv", DataFrame)

---

引入噪声与现实实验 🔬

现在我们已经有了整洁的数据，本节中我们来看看如何模拟现实世界的实验。在现实中，测量数据总是包含噪声。我们可以通过在真实数据上添加随机噪声来模拟这种情况。

假设真实的转换公式是：celsius = (5/9) * fahrenheit - (160/9)。我们生成一些华氏度数据 x，并计算对应的精确摄氏度 y。然后，我们创建带噪声的版本 yy：

using Random
# 真实参数
true_slope = 5/9
true_intercept = -160/9 # 约等于 -17.777...

# 生成数据
x = rand(-10:100, 10) # 10个华氏度数据点
y_true = true_slope .* x .+ true_intercept

# 添加噪声
noise_level = 0.5
yy = y_true .+ noise_level .* randn(10) # 添加正态分布噪声

每次运行实验（即生成一组新的噪声数据），我们拟合出的直线都会略有不同。统计学要解决的核心问题就是：基于这有限的、带噪声的数据，我们对拟合出的直线有多大信心？

---

拟合直线：最小二乘法 📉

上一节我们模拟了带噪声的数据，本节中我们来看看如何从这些数据中找到“最佳”直线。最常用的方法是最小二乘法，它寻找一条直线，使得所有数据点到这条直线垂直距离的平方和最小。

这条蓝色最佳拟合线由斜率 m 和截距 b 定义。我们可以通过线性代数公式计算它们：

# 构建设计矩阵 X：第一列为1（对应截距），第二列为x
X = [ones(length(x)) x]
# 使用线性代数求解最小二乘解：θ = (X'X)^{-1} X'y
θ = (X' * X) \ (X' * yy)
b_estimate, m_estimate = θ[1], θ[2]

也可以用更基础的公式计算，结果相同：

n = length(x)
x_mean = mean(x)
y_mean = mean(yy)

# 计算斜率 m
m_estimate = sum((x .- x_mean) .* (yy .- y_mean)) / sum((x .- x_mean).^2)
# 计算截距 b
b_estimate = y_mean - m_estimate * x_mean

此外，我们还可以估计数据的噪声水平（方差 σ²）：

# 计算残差（数据点到拟合线的垂直距离）
residuals = yy .- (m_estimate .* x .+ b_estimate)
# 估计噪声方差，除数是 n-2（自由度）
sigma_squared_estimate = sum(residuals.^2) / (n - 2)

---

理解线性模型与统计假设 🧠

我们一直在进行拟合，本节中我们来看看这背后完整的统计模型。在标准线性回归中，我们假设数据由以下模型生成：

y = m * x + b + ε

其中：

• x 是已知的自变量（如华氏度）。
• m (斜率) 和 b (截距) 是未知的固定参数。
• ε 是随机噪声，我们假设它服从均值为0、方差为 σ² 的正态分布，即 ε ~ N(0, σ²)。

我们通过实验数据估计出 m_estimate、b_estimate 和 sigma_estimate。由于噪声的存在，每次实验得到的估计值都不同。为了理解这些估计值的可靠性，我们需要进行模拟。

---

通过模拟理解估计量的分布 🔄

统计理论告诉我们，在模型的假设下，斜率 m_estimate 和截距 b_estimate 的估计值本身也服从正态分布。我们可以通过Julia快速模拟成千上万次实验来验证这一点。

以下是模拟10万个实验，并观察截距估计值分布的代码：

function run_one_simulation(x, true_slope, true_intercept, true_sigma)
    # 1. 根据真实模型生成带噪声的y
    y_sim = true_slope .* x .+ true_intercept .+ true_sigma .* randn(length(x))
    # 2. 拟合直线
    X = [ones(length(x)) x]
    θ = (X' * X) \ (X' * y_sim)
    b_sim, m_sim = θ[1], θ[2]
    # 3. 估计噪声
    residuals = y_sim .- (m_sim .* x .+ b_sim)
    sigma_sim = sqrt(sum(residuals.^2) / (length(x) - 2))
    return b_sim, m_sim, sigma_sim
end

# 运行10万次模拟
n_simulations = 100_000
results = [run_one_simulation(x, true_slope, true_intercept, noise_level) for _ in 1:n_simulations]

# 提取所有截距估计值
all_intercepts = [r[1] for r in results]
# 计算其均值和标准差
mean_intercept = mean(all_intercepts)
std_intercept = std(all_intercepts)

模拟会发现：

1. all_intercepts 的分布近似正态。
2. 其均值 mean_intercept 非常接近真实截距 true_intercept。
3. 其标准差 std_intercept 与理论公式计算出的标准误差相符。

对斜率 m_estimate 的模拟会得到类似结论。而对噪声估计 sigma_estimate 的模拟则会发现它服从卡方分布。

---

解读统计输出表格 📊

现在我们已经通过模拟理解了估计量的行为，本节中我们终于可以解读那个“神秘的”统计输出表格了。表格通常包含以下几列：系数、标准误差、t 值和 p 值。

假设我们得到如下表格（以截距为例）：

项	系数估计	标准误差	t 值	p 值
截距	-18.32	0.49	-37.4	<0.001
斜率	0.534	0.033	16.2	<0.001

1. 系数估计： 这就是我们通过最小二乘法计算出的 b_estimate 和 m_estimate。

2. 标准误差： 这是系数估计值的估计标准差。它衡量了如果我们重复实验，系数估计值可能的变化范围。计算公式使用了我们估计的噪声方差 sigma_squared_estimate。
* 截距标准误差 ≈ sigma_estimate * sqrt(1/n + mean(x)^2 / sum((x .- mean(x)).^2))
* 斜率标准误差 ≈ sigma_estimate / sqrt(sum((x .- mean(x)).^2))

3. t 值： 这是最简单的部分，它就是系数估计值除以它的标准误差。
* t = 系数估计 / 标准误差
* 它衡量了系数估计值相对于其波动性的大小。

4. p 值： 它基于 t 分布 计算。我们首先有一个零假设：真实的系数为0（例如，斜率=0意味着x和y没有线性关系）。p值表示，如果零假设成立（系数真为0），我们观察到当前这么大（或更大）的|t|值的概率。
* p值很小（如<0.05）：意味着如果系数真为0，观察到当前数据的情况非常不可能。因此，我们拒绝零假设，认为该系数是“显著的”（很可能不为0）。
* p值较大：则没有足够证据拒绝零假设，该系数可能不显著。

自由度： 在上面的噪声估计公式 sigma_squared_estimate = sum(residuals.^2) / (n - 2) 中，分母 n-2 就是自由度。这里 -2 是因为我们估计了两个参数（斜率和截距），消耗了2个自由度。n-2 也是t分布的自由度参数。

---

总结

本节课中我们一起学习了线性模型与统计模拟的核心内容。

1. 我们首先使用数据框来组织带标签的数据。
2. 通过添加噪声模拟了现实世界的实验，认识到从有限噪声数据中得出的结论具有不确定性。
3. 我们使用最小二乘法拟合直线，得到斜率、截距和噪声的估计值。
4. 我们明确了标准线性模型的统计假设：数据由线性关系加上正态噪声生成。
5. 通过进行大量计算机模拟，我们直观地验证了系数估计量的抽样分布（正态分布），并理解了其波动性。
6. 最后，我们完整解读了统计输出表格：系数是点估计，标准误差衡量其精度，t值是系数与标准误差的比值，而p值则用于检验系数是否显著不为零。

理解这些概念的关键在于，统计学并非提供确凿无疑的答案，而是基于数据和模型假设，对未知世界进行量化的、带有置信度的推断。

L16：优化 🎯

在本节课中，我们将学习优化的基本概念，并通过一个简单的线性拟合问题，探索在Julia中实现优化的多种方法。我们将从精确的数学公式开始，逐步深入到使用优化软件包、梯度计算以及随机梯度下降等现代机器学习中常用的技术。

---

概述

优化是寻找最佳解决方案的过程，在数据科学和机器学习中至关重要。本节课我们将以最小二乘法拟合直线为例，演示如何用多种方式解决同一个优化问题，并比较它们的优缺点。

---

数学定义与精确解

首先，我们明确要解决的问题：给定一组数据点 (x_i, y_i)，找到最佳拟合直线 y = b + m*x，使得所有点的预测值与实际值之差的平方和最小。这个和被称为损失函数。

损失函数的数学公式为：

L(b, m) = Σ (b + m*x_i - y_i)^2

对于这个简单问题，存在精确的解析解。我们可以通过线性代数公式直接计算出最优的截距 b 和斜率 m。

线性代数解法的核心是构建一个设计矩阵 A，其中第一列全为1，第二列为 x 坐标。然后通过求解正规方程 (A^T A) \ (A^T y) 得到参数。

在Julia中，我们可以使用命名元组来清晰地存储和展示结果：

result = (intercept = b, slope = m)

---

使用优化软件包

上一节我们介绍了问题的数学定义和精确解。本节中我们来看看如何使用通用的优化软件包来解决同一个问题。当问题变得复杂，没有简单解析解时，这种方法显得尤为重要。

我们将介绍两个Julia包：Optim.jl 和 JuMP.jl。

使用 Optim.jl

Optim.jl 是一个完全用Julia编写的优化包。我们只需要定义损失函数，并提供一个初始猜测值，它就能自动寻找最小值点。

以下是使用 Optim.jl 进行优化的步骤：

1. 定义损失函数 loss(params)，其中 params 是一个包含 b 和 m 的向量。
2. 调用 optimize(loss, initial_guess)。
3. 从结果中提取优化后的参数。

软件会自动选择算法（如Nelder-Mead）并报告收敛信息，例如函数调用次数和最终精度。

使用 JuMP.jl

JuMP.jl 是一个数学建模语言，它允许用户以更接近数学表达式的形式描述优化问题。它本身不包含求解器，但可以连接多种后端求解器（如Ipopt）。

以下是使用 JuMP.jl 建模的步骤：

1. 创建一个模型 model = Model(Ipopt.Optimizer)。
2. 使用 @variable 定义优化变量 b 和 m。
3. 使用 @objective 定义最小化目标，即我们的损失函数。
4. 调用 optimize!(model) 进行求解。

这种方法对于表达复杂的约束优化问题非常强大。

---

梯度的作用与计算

上一节我们使用了现成的优化包。本节中我们来看看优化背后的一个核心数学概念：梯度。理解梯度有助于我们理解优化器是如何工作的，并让我们能实现自己的简单优化算法。

梯度是一个向量，其每个分量是损失函数对相应参数的偏导数。对于我们的损失函数，梯度指出了参数空间中使函数值上升最快的方向。因此，负梯度方向就是函数值下降最快的方向。

梯度公式如下：

∇L(b, m) = [ ∂L/∂b, ∂L/∂m ] = [ 2*Σ (b + m*x_i - y_i), 2*Σ (b + m*x_i - y_i)*x_i ]

在Julia中，我们有多种计算梯度的方法：

1. 手动推导与编码：根据上面的公式直接编写代码。
2. 有限差分法：通过给参数一个微小的扰动 ϵ，用 (L(param+ϵ) - L(param)) / ϵ 来近似导数。这种方法简单但可能有数值误差，且需要谨慎选择 ϵ。
3. 自动微分：这是现代计算框架的核心技术。它利用链式法则，通过代码本身自动、精确地计算导数，没有符号计算的复杂性和有限差分的误差。在Julia中，可以使用 ForwardDiff.jl 等包轻松实现：

   using ForwardDiff
   gradient = ForwardDiff.gradient(loss, [b, m])
   

自动微分使得“为任意代码计算梯度”成为可能，极大地推动了机器学习和优化领域的发展。

---

实现梯度下降算法

既然我们已经能够计算梯度，本节中我们就可以尝试实现最基本的优化算法——梯度下降。其思想非常直观：就像滑雪者沿着最陡的方向下山一样，我们沿着负梯度方向更新参数。

梯度下降的迭代公式为：

params_new = params_old - η * ∇L(params_old)

其中 η 称为学习率或步长。

以下是实现梯度下降的关键步骤：

1. 初始化参数（例如，设 b=0, m=0）。
2. 循环多次：
   a. 计算当前参数下的梯度 g。
   b. 更新参数：params -= η * g。
3. 循环结束后，params 即为近似的最优解。

算法的成败很大程度上取决于学习率 η 的选择：

• η 太小：收敛速度极慢，需要非常多步迭代。
• η 太大：更新步伐过大，可能无法收敛，甚至在最优点附近震荡或发散。

对于我们的简单问题，我们甚至可以进行精确线搜索，即直接计算出沿梯度方向最优的步长。但在复杂问题中，这通常不可行，因此选择合适的 η 更像一门艺术。

---

随机梯度下降

上一节我们介绍了标准的梯度下降，它需要在每一步计算所有数据点的梯度。本节中我们来看看当数据量非常庞大时，一种更实用的变体：随机梯度下降。

在机器学习中，数据集往往包含数百万甚至数十亿个样本。计算完整梯度（称为“批量梯度下降”）的成本极高。SGD的核心思想是：在每一步迭代中，随机选取一个（或一小批）数据点，只用这个子集的数据来计算梯度并更新参数。

更新公式与梯度下降类似，但梯度 g 是基于随机样本 i 计算的：

∇L_i(b, m) = [ 2*(b + m*x_i - y_i), 2*(b + m*x_i - y_i)*x_i ]
params_new = params_old - η * ∇L_i(params_old)

以下是SGD的特点：

• 优点：每次迭代计算量小，速度快，可以处理海量数据。即使没有收敛到精确解，通常也能快速得到一个足够好的近似解，这在机器学习中往往是可接受的。
• 缺点：更新方向基于噪声估计，路径曲折，需要更多迭代步数才能稳定。对学习率 η 的调整更为敏感。
• 实践：在实际中，我们通常会遍历整个数据集多次（每个循环称为一个“epoch”），并可能随着时间推移逐渐减小 η。

我们通过代码演示了SGD，即使需要上千万次迭代，它最终也能逼近最优解，这体现了其在处理大规模问题时的可行性。

---

高级优化器与总结

在最后一部分，我们回到 Optim.jl 包，看看提供梯度信息如何帮助高级优化器更高效地工作。

当我们能为优化器提供梯度函数时，它可以使用更强大的算法，例如 BFGS 算法。BFGS是一种拟牛顿法，它会构建损失函数曲率的近似，从而智能地调整每一步的更新方向和步长。

在 Optim.jl 中，我们可以这样使用：

result = optimize(loss, gradient, initial_guess, BFGS())

与不提供梯度的Nelder-Mead方法相比，BFGS通常能以少得多的函数调用次数达到相同的精度。这展示了利用问题的一阶信息（梯度）所能带来的巨大效率提升。

---

本节课总结

本节课我们一起学习了优化的基本概念，并以线性拟合为例探索了多种解决途径：

1. 精确解：对于简单问题，直接使用数学公式最快、最准确。
2. 通用优化包：如 Optim.jl 和 JuMP.jl，它们将复杂的优化算法封装成简单接口，适用于没有解析解的问题。
3. 梯度与自动微分：梯度是优化的核心。自动微分技术让我们能轻松为复杂代码计算精确梯度，是当代机器学习框架的基石。
4. 梯度下降：最基本的迭代优化算法，其性能高度依赖于学习率的选择。
5. 随机梯度下降：为大规模数据集设计的优化算法，每次更新只使用一部分数据，牺牲了部分精度以换取计算效率。
6. 高级优化器：如BFGS，在获得梯度信息后，能通过更复杂的数学方法实现快速收敛。

通过这个简单的例子，我们看到了从具体数学解法到通用数值优化，再到适应大数据范式的随机算法的演进过程，这也是现代计算科学解决实际问题的一个缩影。

L17：时间步进与微分方程

在本节课中，我们将学习如何从直观的离散时间模型出发，推导出微分方程，并了解如何使用计算机通过时间步进法来求解微分方程。我们将从简单的灯泡失效模型开始，逐步过渡到更复杂的SIR传染病模型，并探讨数值求解的基本思想。

从离散模型到连续模型

上一节我们介绍了课程将进入新的模块，重点是利用计算思维理解气候科学等领域的模型。这些模型通常由微分方程描述。本节中，我们来看看如何从更直观的离散时间模型自然地引出微分方程。

灯泡失效模型

我们从一个简单的例子开始：建模灯泡（或其他组件）的失效过程。假设我们每天检查一次，记录仍在工作的灯泡数量 N_k，其中 k 代表第几天。每个灯泡每天有概率 p 失效。

以下是模型的核心递推关系：

N_{k+1} = N_k - p * N_k

这表示下一天的数量等于当天的数量减去当天失效的数量。我们可以将其重写，突出变化量：

N_{k+1} - N_k = -p * N_k

这个离散模型可以解析求解，其解为几何衰减：

N_k = N_0 * (1 - p)^k

增加检查频率

现在，我们考虑更频繁地检查，比如每天检查 n 次。假设在每一个长度为 1/n 天的时间段内，失效概率为 p/n。

此时，递推关系变为：

N_{k + 1/n} = N_k * (1 - p/n)

经过一天（即 n 个步骤）后，总变化为：

N_{k+1} = N_k * (1 - p/n)^n

当我们用代码绘制并逐渐增加 n（即减小时间步长）时，会发现离散点构成的曲线逐渐趋近于一条光滑的极限曲线。这暗示着存在一个与离散化细节无关的连续时间模型。

走向连续极限

令时间步长 Δt = 1/n。在步长 Δt 内的递推公式为：

N(t + Δt) - N(t) = -p * Δt * N(t)

两边同时除以 Δt，得到：

( N(t + Δt) - N(t) ) / Δt = -p * N(t)

当 Δt 趋近于0时，左边的表达式就变成了导数 dN/dt 的定义。于是，我们得到了一个微分方程：

dN/dt = -p * N(t)

其初始条件为 N(0) = N_0。

这个微分方程的解是指数衰减函数，验证了离散模型在连续极限下的行为：

N(t) = N_0 * exp(-p * t)

更复杂的模型：SIR传染病模型

上一节我们通过灯泡模型看到了离散与连续的联系。本节中我们来看看一个更复杂、无法解析求解的模型——SIR模型，它常用于描述流行病或谣言传播。

在SIR模型中，人群被分为三类：

• 易感者 (S)：可能被感染的人。
• 感染者 (I)：已患病并可传染他人的人。
• 康复者 (R)：已康复且具有免疫力的人。

总人口数 N = S + I + R。为简化，我们常使用比例 s = S/N, i = I/N, r = R/N。

离散时间SIR模型

假设时间步长为一天。模型基于以下假设：

1. 新感染人数与易感者和感染者的接触成正比（质量作用定律）。
2. 感染者每天以固定概率康复。

以下是离散模型的递推公式：

ΔI = β * s_t * i_t
ΔR = γ * i_t
s_{t+1} = s_t - ΔI
i_{t+1} = i_t + ΔI - ΔR
r_{t+1} = r_t + ΔR

其中 β 是感染率参数，γ 是康复率参数。

连续时间SIR模型

遵循从灯泡模型学到的方法，我们将离散时间步长推广到任意步长 Δt，然后取 Δt → 0 的极限，得到一组耦合的微分方程：

ds/dt = -β * s(t) * i(t)
di/dt =  β * s(t) * i(t) - γ * i(t)
dr/dt =                     γ * i(t)

这是一个非线性系统，因为方程中包含 s(t) * i(t) 这样的乘积项。对于此类模型，通常无法找到解析解（即写出 s(t), i(t), r(t) 关于 t 的显式公式）。

用时间步进法求解微分方程

既然许多有趣的微分方程（如SIR模型）无法解析求解，我们如何利用计算机找到它们的解呢？答案是将微分方程“离散化”，也就是使用时间步进法。

欧拉方法

最基本的时间步进法是欧拉方法。其核心思想是逆向运用我们之前推导微分方程的过程。

对于一个一般形式的微分方程：

dx/dt = f(x)

欧拉方法用有限差分来近似导数：

( x(t + h) - x(t) ) / h ≈ f(x(t))

其中 h 是一个小的时间步长。将其重排，就得到了一个递推公式（时间步进公式）：

x(t + h) = x(t) + h * f(x(t))

或者，如果我们用 x_k 表示第 k 步的值：

x_{k+1} = x_k + h * f(x_k)

对于SIR模型，我们可以将 s, i, r 组成一个向量 x = [s, i, r]，然后同样应用欧拉方法进行迭代计算。

方法的局限性与改进

欧拉方法虽然直观，但精度往往不高，有时甚至无法给出定性正确的解（例如在模拟守恒系统如钟摆时，可能导致能量错误增加）。因此，在实际科学计算中，人们会使用更精密的方法，如龙格-库塔法等。

Julia生态系统提供了一个强大且先进的微分方程求解库 DifferentialEquations.jl（现已整合到 SciML 框架中），其中实现了多种高精度、高效率的数值算法，适用于解决各类复杂的微分方程问题。

总结

本节课中我们一起学习了微分方程与计算思维的核心联系。

1. 我们从离散时间模型（如每天检查灯泡）出发，通过增加检查频率、取连续极限，自然地导出了微分方程。
2. 我们探讨了更复杂的SIR模型，它由一组非线性的耦合微分方程描述，通常无法解析求解。
3. 为了用计算机求解此类方程，我们引入了时间步进法，最基本的是欧拉方法，它将微分方程重新转化为离散递推问题进行近似求解。
4. 我们认识到欧拉方法有其局限性，而在实际应用中，会采用更先进的数值方法，并可以利用像Julia的 DifferentialEquations.jl 这样的专业工具库。

通过理解从离散到连续、再从连续（微分方程）回到离散（数值求解）的这个循环，我们掌握了用计算思维建模和求解连续动态系统的基本流程。

L18：工具库和参数类型 🧰

在本节课中，我们将学习如何使用Julia中的工具库（包）来求解常微分方程，并深入理解这些库是如何通过类型和参数化类型来构建的，以实现高效和灵活的代码。

---

上一节我们介绍了常微分方程的基本概念。本节中，我们来看看如何使用一个专门的Julia包来求解它们，并探索其背后的软件设计原理。

使用微分方程求解包

我们首先需要将微分方程写成一个特定的函数形式。对于一个方程 u̇ = -p * u，其右侧需要定义为一个接受三个参数的函数：变量 u、参数 p 和时间 t。

代码示例：定义微分方程

f(u, p, t) = -p * u

即使方程本身不显式依赖于时间 t，函数定义中仍需包含 t 参数，因为软件包要求函数具有固定的参数结构。

接下来，我们需要指定求解所需的数据：时间区间和初始条件。

代码示例：设置问题参数

tspan = (0.0, 10.0)  # 时间区间
u0 = 100.0           # 初始条件
p = 0.2              # 方程参数

构建与求解ODE问题

Julia的 DifferentialEquations.jl 包围绕类型进行构建。我们使用 ODEProblem 类型来封装微分方程的所有信息。

代码示例：定义并求解问题

using DifferentialEquations
prob = ODEProblem(f, u0, tspan, p)  # 构建问题对象
sol = solve(prob)                    # 求解问题

solve 函数接收这个 ODEProblem 对象，并返回一个 ODESolution 类型的对象，其中包含了求解结果。

处理与可视化结果

solution 对象不仅存储了离散时间点上的解，还包含了进行插值以获取连续解的信息。我们可以像调用函数一样调用它来获取任意时间点的值。

代码示例：使用解对象

sol(3.5)          # 在 t=3.5 处插值求解
plot(sol)         # 绘制解的连续曲线
sol.t             # 获取求解所用的时间点向量
sol.u             # 获取对应时间点的解值向量

plot(sol) 能够绘制出平滑曲线，这得益于包中定义的绘图配方，它告诉绘图库如何绘制这种特定类型的对象。

求解方程组

该方法同样适用于方程组。例如，对于SIR传染病模型，我们将三个变量放入一个向量中，并定义一个返回向量值的函数。

代码示例：定义SIR模型方程组

function sir!(du, u, p, t)
    s, i, r = u
    β, γ = p
    du[1] = -β * s * i
    du[2] = β * s * i - γ * i
    du[3] = γ * i
end

u0 = [0.99, 0.01, 0.0]  # 初始向量 [S, I, R]
p = [0.3, 0.1]          # 参数 [β, γ]
tspan = (0.0, 50.0)
prob = ODEProblem(sir!, u0, tspan, p)
sol = solve(prob)
plot(sol)

探索类型与参数化类型

包中的复杂功能是通过精心设计的类型系统实现的。例如，ODEProblem 类型内部包含多个字段，其具体类型由参数化类型决定，这确保了代码的效率和专一性。

以下是参数化类型的一个简单示例：

代码示例：自定义参数化类型

struct MyODEProblem{T, S}
    t0::T
    tfinal::T
    u0::S
end

# 创建不同类型的实例
prob1 = MyODEProblem(3, 4, 100.0)     # 类型为 MyODEProblem{Int64, Float64}
prob2 = MyODEProblem(3.0, 4.0, 100.0) # 类型为 MyODEProblem{Float64, Float64}
prob3 = MyODEProblem(0.0, 10.0, [1.0, 2.0]) # 类型为 MyODEProblem{Float64, Vector{Float64}}

通过参数 T 和 S，同一个结构体可以灵活地处理不同类型的数据（如标量或向量），同时Julia编译器能为每种类型组合生成高效的特化代码。

可调用对象

我们可以让自定义类型的对象像函数一样被调用。这是通过为类型定义函数调用行为实现的。

代码示例：定义可调用对象

struct SimpleEulerOutput{T, S}
    times::Vector{T}
    values::Vector{S}
end

# 定义调用行为：线性插值
function (sol::SimpleEulerOutput)(t)
    # 在此实现线性插值逻辑，使用 sol.times 和 sol.values
    return "在时间 $t 处插值"
end

sol_obj = SimpleEulerOutput([1, 2], [3, 4])
sol_obj(3.5)  # 输出：在时间 3.5 处插值

---

本节课中我们一起学习了如何使用 DifferentialEquations.jl 包来求解常微分方程，理解了该包如何利用类型系统（如 ODEProblem 和 ODESolution）来组织代码和数据。我们还探讨了参数化类型如何提供灵活性并保证性能，以及如何创建可调用对象。这种基于类型的设计模式是许多高性能Julia包的核心，它连接了数学问题表述与高效的软件实现。

L19：为什么我们不能预测天气 🌦️

在本节课中，我们将探讨一个核心问题：为什么我们无法精确预测天气，却可以对气候做出预测。我们将通过建立简单的数学模型，理解天气系统的复杂性和混沌本质，并学习动力系统中的关键概念，如稳定性、分岔和混沌。

天气与气候的区别

上一节我们介绍了课程主题，本节中我们来看看天气与气候的基本区别。

天气指的是短期内的大气状态，例如温度、降水、风速等。我们可以通过观测当前云层和气压，并利用描述空气和水运动规律的偏微分方程模型，来预测未来几小时甚至几天的天气。这些模型基于我们已充分理解的物理过程。

气候则是指长期的平均天气状况，例如一个地区的年平均温度。虽然每天的天气变化莫测，但当我们对时间（如整个三月份）和空间进行平均后，气候模式就变得相对稳定和可预测。气候变化指的是这种长期平均状态的缓慢改变。

从简单模型开始：逻辑斯蒂方程

为了理解复杂系统，我们从一个简单的模型开始。这个模型描述细菌在有限资源下的种群增长。

我们用一个变量 x 表示细菌数量（相对于环境承载能力的比例）。其增长规律可以用以下常微分方程描述：
dx/dt = x * (1 - x)
其中，x * (1) 项代表指数增长趋势，- x^2 项代表资源有限导致的竞争和饱和效应。

以下是使用数值方法（如欧拉法）模拟该方程行为的步骤：

1. 定义微分方程的右侧函数 f(x) = x * (1 - x)。
2. 选择一个初始值 x0、时间步长 h 和总模拟时间 T。
3. 迭代计算：x_{n+1} = x_n + h * f(x_n)。

模拟结果显示：

• 如果初始值 x0 在 0 和 1 之间，种群数量会增长并最终稳定在 x = 1（即达到环境最大承载量）。
• 如果初始值 x0 = 1，种群数量保持不变。
• 如果初始值 x0 = 0，种群数量也保持不变。
• 如果初始值 x0 > 1，种群数量会下降至 1。
• 如果初始值 x0 < 0（虽无物理意义），种群数量会趋向负无穷。

x = 0 和 x = 1 都是系统的不动点（Fixed Points），即状态不再随时间变化的点。x = 1 是稳定不动点，因为附近的轨迹都会被吸引过来；x = 0 是不稳定不动点，因为附近的轨迹会远离它。

我们可以不通过模拟，而通过绘制相图（Phase Portrait）来定性分析系统行为。相图是在状态空间（这里是 x 轴）上，用箭头表示每个点处 dx/dt 的方向（即变化趋势）。箭头指向 x = 1，说明它是吸引子；箭头远离 x = 0，说明它是排斥子。

引入参数：分岔现象

上一节我们看了一个固定方程的行为，本节中我们来看看当方程中包含可调参数时，系统行为如何发生质变。

考虑方程：dx/dt = μ + x^2
这里 μ 是一个参数。不动点满足 μ + x^2 = 0。

当我们改变 μ 时：

• 当 μ < 0 时，方程有两个不动点：一个稳定，一个不稳定。
• 当 μ 增大时，两个不动点逐渐靠近。
• 当 μ = 0 时，两个不动点碰撞并消失。
• 当 μ > 0 时，方程没有实不动点，所有轨迹都趋向无穷。

这种因参数微小连续变化而导致系统定性行为（如不动点数量或稳定性）发生突然改变的现象，称为分岔（Bifurcation）。μ = 0 是分岔点。

将不同 μ 值下的相图堆叠起来，就得到了分岔图，它能清晰展示系统行为随参数的变化。

更复杂的动态：双稳性与滞后

现在考虑一个三次方程：dx/dt = μ * x - x^3
其不动点满足 x * (μ - x^2) = 0。

分岔图显示：

• 当 μ ≤ 0 时，只有一个稳定不动点 x = 0。
• 当 μ > 0 时，有三个不动点：x = 0（不稳定），x = ±√μ（稳定）。

在 μ > 0 的区域，系统具有双稳性（Bistability），即存在两个稳定的不动点。系统最终到达哪一个，取决于初始条件。

如果缓慢地改变参数 μ：

• 从 μ < 0 开始增加，系统状态会沿着 x = 0 的稳定分支移动。
• 当 μ 越过 0 进入正区域后，x = 0 变得不稳定，系统状态会突然跳跃到 x = √μ 的稳定分支上。
• 如果此后减小 μ，系统状态会沿着 x = √μ 分支移动，而不会在 μ 刚小于 0 时就跳回 x = 0 分支，直到 μ 进一步减小到另一个临界值才会跳回。

这种系统状态变化路径依赖于历史（参数变化方向）的现象，称为滞后（Hysteresis）。某些简单气候模型中就存在这种双稳性和滞后现象。

二维系统：振荡与极限环

在一维系统中，轨迹只能趋向或远离不动点。在二维或更高维系统中，会出现更丰富的行为，例如周期性振荡。

以布鲁塞尔子模型（Brusselator）为例，这是一个描述化学振荡反应（如BZ反应）的模型：
dx/dt = a - (b+1)*x + x^2 * y
dy/dt = b * x - x^2 * y
其中 x, y 是化学物质浓度，a, b 是参数。

在二维相平面 (x, y) 上模拟该系统：

• 对于某些参数，所有轨迹都螺旋式收敛到一个稳定的不动点（阻尼振荡）。
• 当参数 b 超过某个临界值（发生霍普夫分岔，Hopf Bifurcation）时，不动点失去稳定性，系统产生一个稳定的极限环（Limit Cycle）。
• 极限环是一个孤立的闭合轨迹，附近的轨迹都会被吸引到它上面，系统进入稳定的周期性振荡状态。这是二维系统特有的行为。

三维系统：混沌与洛伦兹吸引子

最后，我们进入模拟天气的核心——洛伦兹系统。它是由气象学家爱德华·洛伦兹从简化的大气对流模型中推导出的三维常微分方程组：
dx/dt = σ * (y - x)
dy/dt = x * (ρ - z) - y
dz/dt = x * y - β * z
其中 σ, ρ, β 是参数，通常取 σ=10, β=8/3，而 ρ 作为关键参数变化。

当 ρ 较小时，系统行为简单，轨迹趋向于一个或两个不动点。随着 ρ 增大，系统经历一系列分岔。当 ρ 达到约 28 时，系统出现令人惊奇的行为：

• 轨迹不再稳定于不动点或周期轨道，而是在一个复杂的、类似蝴蝶翅膀形状的集合上无限游走，这个集合称为洛伦兹吸引子。
• 该吸引子具有分形结构，且系统表现出确定性混沌。

混沌的核心特征是对初始条件的极端敏感性：两个无限接近的初始点，其轨迹会随时间呈指数级分离。我们可以计算两条初始距离为 ε 的轨迹之间的距离 d(t)，会发现 d(t) ~ ε * e^(λt)，其中 λ > 0 称为李雅普诺夫指数。这意味着即使模型是完全确定的，长期的精确预测也变得不可能，因为初始测量任何微小的误差都会被迅速放大。这从根本上解释了为何天气预报在超过一两周后就会失效。

然而，虽然单个轨迹不可长期预测，但混沌系统在吸引子上的长期统计平均性质（如平均值、方差）却是稳定的。对于洛伦兹系统，即使两条轨迹完全不同，它们各自在长时间内的坐标平均值也非常接近。这正对应了“气候”（长期统计平均）的可预测性，尽管“天气”（具体瞬时状态）是混沌的。

总结

本节课中我们一起学习了动力系统的基本概念及其在理解天气与气候预测问题中的应用。

我们首先区分了天气（短期、具体、混沌）和气候（长期、平均、可预测）。然后从简单的一维逻辑斯蒂方程入手，学习了不动点、稳定性、相图等概念。通过引入参数，我们观察了分岔现象，以及双稳性和滞后等非线性效应。

进入二维系统，我们看到了极限环和周期性振荡如何从霍普夫分岔中产生。最后，在三维洛伦兹系统中，我们遇到了确定性混沌、对初始条件的敏感性以及奇怪的吸引子。这揭示了天气不可长期预测的数学根源。同时，我们也认识到，混沌系统的长期统计行为可能是稳定的，这为气候预测提供了可能性。

这些非线性动力学概念不仅适用于气象学，也广泛存在于物理、化学、生物乃至社会科学等众多领域。理解这些思想，是培养计算思维、应对复杂系统挑战的重要一步。

L20：第一个气候模型 🌍

在本节课中，我们将学习如何构建一个简单的气候模型。我们将从基础的微分方程入手，逐步引入太阳辐射、地球热辐射以及温室效应等物理概念，最终整合成一个能够模拟地球温度变化的“零维”能量平衡模型。通过这个过程，你将理解计算思维如何帮助我们直观地探索和理解复杂的数学模型。

---

第一部分：从微分方程开始 📈

上一节我们介绍了课程的主题。本节中，我们来看看如何用计算思维来理解微分方程，这是构建气候模型的基础。

我们从一个关于时间 t 的函数 y(t) 开始。设初始时间 t = 0，初始值为 y0。我们用 y‘ 表示 dy/dt。

在传统的微分方程课程中，我们可能会从简单的方程开始，例如 y‘ = a（a 为常数），并求得其符号解 y(t) = a*t + y0。

接下来，我们增加一个线性项，考虑方程 y‘ = a - b*y。这个方程有精确的符号解。当 t 趋于无穷大时，y(t) 会趋近于一个平衡点 a/b。这个平衡点可以通过令 y‘ = 0 求得。

我们还可以加入一个与时间相关的强迫项 f(t)，方程变为 y‘ = a - b*y + f(t)。其解包含一个积分项。虽然这些符号解包含了所有信息，但它们通常不够直观。

我更喜欢通过数值求解和可视化来理解微分方程。Julia 语言拥有由 Christopher Rackauckas 编写的卓越的 ODE 求解器套件。使用方法如下：

1. 定义函数 f(y, p, t)，描述 y‘。
2. 使用 ODEProblem 设置问题，包括初始条件、时间范围和参数。
3. 使用 solve 命令求解。
4. 对结果进行绘图和分析。

通过交互式地调整参数（如 a, b, c），并观察方向场和解曲线的变化，我们可以更直观地理解每个项在方程中的作用：

• 参数 a：作为常数项，它均匀地影响所有位置，推动解曲线整体上升。
• 参数 b：与 y 线性相关（-b*y），它促使解曲线趋向于一个平衡点，防止其无限增长。
• 强迫项 f(t)：作为时间的函数，它会改变平衡点的位置或使解曲线随时间产生特定变化。

这种通过计算和可视化来探索方程行为的方式，正是计算思维的核心——利用计算机作为理解数学的工具。

---

第二部分：构建零维气候模型 ☀️❄️

上一节我们介绍了如何用计算思维理解微分方程。本节中，我们来看看如何将这些方程应用到气候科学中，构建一个简单的“零维”能量平衡模型。该模型将地球视为一个点，只考虑全球平均温度。

气候模型主要考虑三个物理过程：

1. 吸收的太阳辐射：来自太阳的能量，使地球变暖。
2. 向外热辐射：地球表面向太空辐射的能量，使地球冷却。
3. 人为温室效应：大气中温室气体（如二氧化碳）捕获部分本应逃逸到太空的热辐射，导致增温。

以下是这三个过程对应的微分方程核心项：

1. 仅考虑吸收的太阳辐射
这对应于最简单的方程：温度变化率 ∝ 常数。
吸收的太阳辐射功率 S 约为 1368 W/m²。地球反射一部分太阳光，反射率 α 约为 0.3。同时，由于地球是球体，接收太阳辐射的有效面积是其横截面积（πR²），而非总面积（4πR²），因此需要除以 4。
此外，海洋的热容量 C（约为 51 J/m²/°C）会减缓温度变化。
因此，温度 T 的变化方程是：
C * dT/dt = S * (1 - α) / 4
代入数值，右端约为 239.4 W/m²。如果地球只吸收太阳热量而不散发，温度将直线上升，在工业革命（1850年）至今的170多年里，会达到无法想象的高温。

2. 加入向外热辐射
这对应于方程 y‘ = a - b*y 的形式。向外热辐射近似为温度的线性函数：b * (T0 - T(t))，其中 T0 是平衡温度，负号表示该效应使温度恢复平衡。
我们利用工业革命前的平均温度 T0 = 14°C 作为平衡温度来确定参数关系。通过拟合数据，可以确定系数 b 约为 1.3 W/m²/°C。
现在方程变为：
C * dT/dt = S*(1-α)/4 - b*(T - T0)
在这个模型中，无论初始温度如何，系统最终都会稳定在 14°C 的平衡状态，这模拟了地球自然的温度调节能力。

3. 加入人为温室效应
这对应于增加一个强迫项 f(t)。温室效应强度与大气中二氧化碳（CO₂）浓度的对数成正比。
公式为：F * log₂( CO₂(t) / CO₂_pre )，其中 F 是强迫系数（约 5 W/m²），CO₂_pre 是工业革命前的 CO₂ 浓度（约 280 ppm）。
因此，完整的零维模型方程为：
C * dT/dt = S*(1-α)/4 - b*(T - T0) + F * log₂( CO₂(t) / CO₂_pre )
现在，温度的变化不仅受自然过程影响，也受人类活动导致的 CO₂ 浓度变化 CO₂(t) 驱动。

---

第三部分：使用真实数据与预测未来 📊🔮

上一节我们建立了气候模型的微分方程。本节中，我们来看看如何利用真实数据校准模型，并用它来预测未来的可能情景。

首先，我们需要知道历史 CO₂ 浓度数据 CO₂(t)。我们可以从美国夏威夷莫纳罗亚观测站获取自1958年以来的月度数据。在 Julia 中，我们可以使用 CSV.read 和 DataFrame 来加载和处理这些数据。数据中可能包含缺失值（如 -99.99），我们可以通过布尔掩码进行过滤，只选取有效数据。

# 示例：使用掩码过滤有效数据
valid_mask = data[!, :co2] .> 0
valid_dates = data[valid_mask, :date]
valid_co2 = data[valid_mask, :co2]

将有效数据绘制出来，可以看到著名的“基林曲线”，它清晰地显示了大气 CO₂ 浓度的持续上升趋势（从约315 ppm升至现今超过410 ppm）。我们可以用一个多项式（例如三次函数）来拟合这条曲线，得到一个平滑的 CO₂(t) 函数用于模型。

其次，我们有全球平均温度的历史观测数据（例如来自 NASA）。将我们模型计算出的温度曲线与真实观测数据进行比较，可以评估模型的准确性。通过调整参数 b（热辐射反馈强度）和 C（海洋热容量），我们可以获得对历史数据更好的拟合。这个过程通常涉及最小二乘法等优化技术。

最后，基于不同的未来 CO₂ 排放情景，我们可以运行模型进行预测。例如：

• 低排放情景：假设人类采取有力措施，将 CO₂ 浓度控制在500 ppm以下。
• 高排放情景：假设排放不受控制，CO₂ 浓度可能升至1200 ppm。

将这两种情景的 CO₂(t) 函数代入我们的微分方程并求解，就可以得到对应的未来温度预测轨迹。《巴黎协定》的目标是将全球升温控制在比前工业化时期高2°C以内。模型预测显示，低排放情景有望接近这一目标，而高排放情景将导致远超2°C的升温，带来严重后果。

---

总结

本节课中，我们一起学习了如何从简单的微分方程出发，逐步构建一个“零维”地球能量平衡气候模型。我们首先通过计算思维和可视化工具直观理解了微分方程中各项参数的作用。然后，我们将吸收太阳辐射、向外热辐射和人为温室效应这三个物理过程转化为数学模型中的具体项。最后，我们利用真实的 CO₂ 和温度观测数据来校准模型，并基于不同的排放情景对未来温度变化进行了预测。

这个简单模型虽然忽略了地理差异、云层、海洋环流等复杂细节（这些需要高性能计算机运行复杂模型），但它清晰地揭示了气候变化背后的核心机制：人类活动导致的温室气体增加，正在破坏地球原有的热平衡。理解这个基本框架，是每个人认识气候变化科学起点的重要一步。

L21：如何在软件上进行协作 👥

在本节课中，我们将学习如何在软件项目中进行协作。我们将重点介绍一个名为GitHub的工具，它被全球开发者广泛用于代码协作。我们将了解它与日常工具（如Google Drive）的区别，并亲自动手创建一个GitHub仓库、提交更改，甚至为一个开源项目做出贡献。

概述：为何需要专门的协作工具？

在开始之前，我们先思考一个简单的问题：最简单的在线协作方式是什么？答案是互相发送电子邮件附件。这对于小型项目或学校论文来说可能有效，但对于软件项目，我们需要更强大的工具。软件项目有其特殊性：微小的代码改动可能导致整个程序崩溃；多人需要同时在不同分支上工作；开发者热衷于自动化流程。这些需求催生了像GitHub这样的版本控制系统。

第一部分：为何不使用Google Drive进行软件协作？ 🤔

上一节我们提到了软件协作的特殊性，本节中我们来看看为什么像Google Drive这样的实时同步工具不适合软件开发。主要有三个原因：

1. 细粒度同步导致程序频繁崩溃：在软件中，一个字符的错误就可能导致整个项目无法运行。如果像Google Drive那样同步每一次按键，那么在修改代码的中间状态（例如，将 sqrt 改为 log 的过程中），项目会多次处于“损坏”状态。我们更希望手动控制何时将一组完整的、可工作的更改发布出去。
2. 需要分支与合并功能：在软件开发中，经常需要创建代码的独立副本（称为“分支”）来尝试新功能或修复错误，而不会影响主版本。之后，需要将不同分支上的更改合并到一起。Google Drive的文件副本机制很难优雅地处理这种复杂的合并操作。
3. 热爱自动化：开发者喜欢自动化一切。GitHub集成了强大的自动化工作流（如自动测试、自动部署网站），这远超出了普通文档协作工具的范围。这种自动化（常被称为DevOps）是软件开发效率的关键，但也增加了复杂性。

第二部分：为开源项目做出贡献 🌍

了解了基本概念后，我们来看看如何实际参与一个开源项目。本节我们将学习如何通过提交“拉取请求”来改进一个气候模型的文档。

以下是向开源项目提交文档更改的基本步骤：

1. 找到目标：浏览项目文档或代码，发现可以改进的地方，例如一个拼写错误或表述不清的句子。
2. 定位代码：在项目的GitHub仓库中找到定义该文档的源文件。你可以使用仓库的搜索功能来查找特定文本。
3. 编辑文件：在GitHub网站上，点击文件右上角的编辑按钮（铅笔图标）。系统可能会提示你需要在一个分支上操作，通常选择默认分支即可。
4. 提交更改：修改内容后，在页面底部填写更改描述，并选择“创建新分支并提交拉取请求”。
5. 发起拉取请求：提交后，会自动创建一个“拉取请求”。项目维护者会审查你的更改，可能会提出修改意见。你可以根据反馈继续更新这个拉取请求中的代码，直到它被合并到主项目中。

核心概念：拉取请求 是你向项目维护者提出的、希望将你的代码更改合并到主代码库的正式请求。

第三部分：创建并管理你自己的GitHub仓库 🛠️

现在，让我们创建自己的代码仓库，体验从本地到远程的完整协作流程。

1. 创建仓库：登录GitHub，点击右上角“+”号，选择“New repository”。为仓库命名（例如 hello-mit），选择公开或私有，建议初始化一个README文件。
2. 克隆到本地：为了在本地电脑上工作，你需要“克隆”这个仓库。推荐使用 GitHub Desktop 这类图形化工具。在仓库页面点击“Code”按钮，选择“Open with GitHub Desktop”，并指定本地存放路径。
3. 进行本地更改：用你喜欢的编辑器（如VS Code）打开本地仓库文件夹，修改或创建文件。例如，修改 README.md 文件。
4. 提交更改：在GitHub Desktop中，你会看到文件更改的差异。勾选想要提交的文件，编写清晰的提交描述，然后点击“Commit”。
5. 推送到远程：提交只是将更改记录在本地。需要点击“Push origin”才能将本地提交同步到GitHub网站上的远程仓库。
6. 获取更新：在协作中，他人也可能修改了远程仓库。经常点击“Fetch origin”来获取远程的最新更改，避免冲突。

第四部分：处理常见问题与冲突 ⚠️

在协作中，当你和他人的修改冲突时，Git会阻止自动合并，要求你手动解决。这对于新手可能有些棘手。

以下是解决冲突的一个简单策略：

1. 备份你的工作：将你修改过的、但尚未成功提交/推送的文件复制到仓库外的安全位置。
2. 重置本地仓库：在GitHub Desktop中，你可以右键点击更改的文件选择“Discard changes”，或者更彻底地删除整个本地仓库文件夹。
3. 重新克隆：从GitHub网站重新克隆一份最新的仓库到本地。
4. 重新应用更改：将之前备份的修改，手动合并到新克隆的文件中，小心处理冲突的部分。
5. 提交并推送：完成合并后，重新提交并推送你的更改。

注意：这是初学者的权宜之计。随着经验增长，你应该学习使用分支、git stash或IDE内置的合并工具来更优雅地解决冲突。

第五部分：测试与更多贡献方式 🧪

为代码编写测试是保证软件质量的重要实践。测试就像一份保险，确保未来的修改不会意外破坏现有功能。

测试驱动开发 是一种先写测试，再写实现代码的方法。这能帮助你明确目标，并产生更健壮的代码。

# 一个简单的测试示例
function double(x)
    return x * 2
end

# 测试用例
@assert double(2) == 4
@assert double(3) == 6
@assert double(0) == 0

除了提交代码，为开源项目做贡献的方式还有很多：

• 报告问题：详细描述你遇到的Bug。
• 改进文档：修正错误或补充示例。
• 提交最小可复现案例：当发现Bug时，提供一个能重现问题的最简代码片段，极大帮助维护者定位问题。

总结

本节课中我们一起学习了软件协作的核心工具GitHub。我们了解了它与普通文档协作工具的区别，动手向开源项目提交了更改，创建并管理了自己的代码仓库，并探讨了处理冲突的基本方法和测试的重要性。记住，开始使用Git时可能会感到复杂，但不要畏惧。从创建一个个人仓库开始，大胆尝试提交、分支等操作，这是掌握这一强大协作工具的最佳途径。

L22：雪球地球与建模 🌍❄️

在本节课中，我们将学习一个重要的数学现象——滞后现象，并探索它与地球气候模型，特别是“雪球地球”假说之间的联系。我们将通过求解微分方程、分析多稳态系统，并最终理解地球气候如何可能在不同状态间切换。

概述：从数学现象到气候模型

上一节我们介绍了气候建模的基本概念。本节中，我们将首先深入探讨一个关键的数学概念——滞后现象，然后将其应用于理解地球气候历史中一个引人入胜的假说：“雪球地球”。

第一部分：数学基础——符号函数与不连续性

为了理解后续的模型，我们首先需要认识一个基础函数：符号函数。

在数学中，符号函数 sign(x)（有时写作 sgn(x)）定义如下：

• 当 x > 0 时，sign(x) = 1
• 当 x = 0 时，sign(x) = 0
• 当 x < 0 时，sign(x) = -1

这是一个在 x=0 处不连续的简单函数。在计算机的IEEE标准中，像 1/0 这样的运算会得到 Inf（正无穷），-1/0 会得到 -Inf（负无穷），它们是有效的数值表示。

第二部分：探索滞后现象

我们将使用计算思维来探索“滞后现象”。在物理学（如磁性研究）中，滞后现象指系统的状态不仅取决于当前条件，还取决于它到达当前状态的历史路径。

为了演示这一点，我们研究一个简单的函数族，它由符号函数构成，并带有一个参数 a：

f(y) = sign(y) + a - y

这个函数图像看起来像两条平行的直线，在 y=0 处有一个跳跃间断。

以下是该函数根（即 f(y)=0 的解）的情况分析：

• 当 a < 1 时，y = a - 1 是一个根（对应负值部分）。
• 当 a > -1 时，y = a + 1 是一个根（对应正值部分）。
• 当 -1 < a < 1 时，函数有两个根。

我们可以将这个函数视为三维空间中的一个曲面（z = f(y, a)），它由两个平行的半平面组成。该曲面与平面 z=0 相交，会产生两条线，这直观地展示了参数 a 如何影响根的数量和位置。

第三部分：微分方程中的多稳态与滞后

现在，我们将 f(y) 作为微分方程 dy/dt = f(y) 的右侧。通过数值求解，我们可以观察系统的行为。

• 当 a 远小于 -1 时，只有一个稳定的平衡点（负值）。
• 当 a 介于 -1 和 1 之间时，存在两个稳定的平衡点（一正一负）和一个不稳定的平衡点。
• 当 a 远大于 1 时，只有一个稳定的平衡点（正值）。

关键现象出现在 -1 < a < 1 的区间。假设系统从 a 值很小（如 -3）开始，并处于对应的负平衡点。然后我们缓慢地增加 a。

1. 系统会连续地沿着负值平衡点的路径移动，即使 a 已经变为正数，它仍然会停留在负值状态。
2. 只有当 a 增加到超过 1 时，负值平衡点消失，系统才会突然跳跃到正值平衡点。
3. 如果我们随后开始缓慢减小 a，系统会沿着正值平衡点的路径移动，直到 a 减小到低于 -1 时，才会再次突然跳回负值状态。

这个循环揭示了滞后现象：对于同一个 a 值（例如 a=0），系统可能处于两个不同的稳定状态（正或负），具体取决于它是从更热还是更冷的历史状态变化而来的。系统的当前状态“记得”它的过去。

第四部分：应用于气候——能量平衡模型

上一讲我们介绍了一个简单的能量平衡模型：地球温度的变化率 = 吸收的太阳辐射 - 向外散发的热辐射。

本节课我们简化模型，忽略人类活动产生的二氧化碳影响，专注于史前气候。关键的变化是：我们不再将地球反照率 α（反射率）视为常数。冰面比水面或陆地反射更多的阳光（反照率更高）。

因此，我们引入一个依赖于温度 T 的反照率函数 α(T)：

• 当温度很高（无冰）时，α 较低（例如 0.3），吸收大部分太阳能。
• 当温度很低（全球冰封）时，α 较高（例如 0.5），反射大部分太阳能。
• 在中间温度，存在一个平滑的过渡区。

将这个 α(T) 代入能量平衡方程后，我们得到了一个与之前数学例子结构相似的方程。绘制“吸收的太阳辐射”和“向外散发的热辐射”随温度变化的曲线，我们会发现它们存在三个交点，即三个平衡点：

1. 一个高温稳定平衡点（约 14°C，对应无冰或少量冰的“间冰期”地球）。
2. 一个低温稳定平衡点（约 -40°C，对应完全冰封的“雪球地球”）。
3. 一个中间的不稳定平衡点（约 -7.5°C）。

第五部分：雪球地球假说与滞后

地质证据表明，地球历史上可能至少发生过三次全球性的冰封事件，即“雪球地球”。我们的模型为这种现象提供了可能的解释。

根据滞后现象的原理：

• 当地球处于温暖的间冰期状态，即使太阳亮度缓慢降低（相当于模型中的 a 减小），气候系统也会倾向于保持在温暖平衡点附近，直到超过某个临界阈值。
• 一旦越过阈值，系统会迅速切换到雪球地球的稳定状态。
• 同样，要逃离雪球地球状态，仅仅依靠太阳亮度缓慢增加是不够的。模型显示，即使将太阳亮度增加到远超今日水平，系统仍可能停留在冰冻状态。这引出了主流理论：是火山持续喷发的大量二氧化碳（温室气体）积累了数百万年，最终足以克服高反照率的冷却效应，将地球推过另一个临界点，使其突然变暖，回到间冰期状态。

更复杂的模型（如考虑纬度分布的“水行星”模型）可以模拟冰盖从极地向赤道扩展的过程，展示了从部分冰封到全球冰封的突变。

总结与启示

本节课我们一起学习了：

1. 滞后现象的数学本质：多稳态系统中，状态取决于历史路径。
2. 如何通过简单的微分方程和计算机模拟来可视化和理解这一现象。
3. 将滞后现象应用于一个简化的地球能量平衡气候模型。
4. 用该模型解释了“雪球地球”假说中气候状态如何发生突变式切换，以及逃离全球冰封可能需要温室气体的长期积累。

核心启示是：地球气候系统并非总是线性、渐变响应外界变化。它可能存在稳定的状态和关键的“引爆点”。当前人类活动在极短时间内大幅增加温室气体浓度，正是在以一种史无前例的速度推动气候系统，这提醒我们不应将当前气候的稳定性视为理所当然。通过计算建模分离和探究这些复杂现象，我们能获得比传统方法更深刻的洞察力。

L23：一维平流扩散偏微分方程 🧮

在本节课中，我们将学习偏微分方程的基础知识，特别是与气候模型相关的平流和扩散过程。我们将从计算思维的角度出发，理解如何通过离散化和数值方法来求解这些方程，并最终将它们结合起来模拟物理现象。

概述 🌍

之前，我们将地球温度建模为一个随时间变化的标量。然而，现实中的温度不仅随时间变化，也随空间位置（如纬度）变化。为了更精确地建模，我们需要引入依赖于空间和时间的温度场 T(x, t)。温度的变化源于热量的流动，这涉及到两种主要物理过程：平流（热量随流体运动）和扩散（热量从高温区域向低温区域自发传播）。本节课的目标是理解并数值求解描述这些过程的偏微分方程。

从离散化开始 📊

为了在计算机上处理连续的温度场，我们首先需要将其离散化。我们将空间区间（例如，从赤道到极点）划分为 N 个等宽的小单元，每个单元的宽度为 Δx。在每个单元的中心点（称为网格点），我们记录该点的温度值 T_i^n，其中 i 是空间索引（表示第 i 个单元），n 是时间步长索引（表示第 n 个时间点）。这样，连续的场 T(x, t) 就被近似为一系列离散值 T_i^n。

平流过程 🌊

平流描述了物理量（如热量或污染物浓度）随流体整体运动而被“携带”的过程。想象一条河流，河水中某处的温度分布会随着水流整体向右移动。

物理图像与离散模型

考虑一个流体单元，其温度为 T_i^n。在时间步长 Δt 内，流体以速度 u 向右运动。因此，一部分原本位于左侧相邻单元 i-1 中的“热量”会流入当前单元 i，同时当前单元中的一部分热量也会流出到右侧单元 i+1。

流入当前单元的热量比例约为 (u Δt) / Δx。基于此质量守恒原理，我们可以写出当前单元在下一时间步的温度更新公式：

T_i^{n+1} = T_i^n + (u Δt / Δx) * (T_{i-1}^n - T_i^n)

这个公式就是平流过程的离散时间步进方案。它表明，下一时刻 i 点的温度等于当前时刻的温度，加上从左边邻居流入的热量，减去从自身流向右边的热量。

边界条件与代码实现

在计算时，对于两端的单元（如 i=1），其左侧邻居 i-1=0 并不存在。我们需要定义边界条件。一种简便的方法是采用周期性边界条件，即将空间域首尾相连，形成一个环。这样，第一个单元的左侧邻居就是最后一个单元。

以下是平流过程的 Julia 代码实现核心部分：

function advection_step(T, u, Δt, Δx)
    N = length(T)
    T_next = similar(T) # 创建新数组存储下一步结果
    α = u * Δt / Δx # 计算 Courant 数

    # 内部单元更新
    for i in 2:N-1
        T_next[i] = T[i] + α * (T[i-1] - T[i])
    end

    # 周期性边界条件处理
    T_next[1] = T[1] + α * (T[N] - T[1])   # 第一个单元的左邻居是最后一个单元
    T_next[N] = T[N] + α * (T[N-1] - T[N]) # 最后一个单元的右邻居是第一个单元（在下一轮计算中体现）

    return T_next
end

从离散到连续：平流方程

如果我们对上述离散方程进行重新排列，并考虑当 Δt 和 Δx 都趋近于零的极限情况，就可以推导出连续的平流方程：

∂T/∂t = -u * ∂T/∂x

这里，∂T/∂t 是温度随时间的变化率，∂T/∂x 是温度在空间上的变化率（梯度）。该方程表明，某点温度随时间的变化，与该点温度在空间上的梯度成正比，方向与流速 u 相反。

数值方法的挑战

直接使用上面推导的简单离散格式（称为迎风格式）进行模拟时，可能会引入非物理的数值扩散，即温度轮廓在平移的同时会变得平滑、振幅衰减。为了更精确地模拟纯平流，通常需要使用更复杂的离散格式，例如中心差分格式：

T_i^{n+1} = T_i^n - (u Δt / (2Δx)) * (T_{i+1}^n - T_{i-1}^n)

这种格式同时考虑了左右邻居的信息，在满足稳定性条件 (|u|Δt/Δx ≤ 1，即 CFL 条件) 时，能更好地保持解的形态。

上一节我们介绍了平流过程，它描述了物理量随流体的整体运动。接下来，我们看看另一种基础物理过程——扩散。

扩散过程 🔥

扩散描述了物理量（如热量、浓度）从高值区域向低值区域自发传播，以达到均匀分布的趋势。即使流体静止，扩散也会发生，例如房间内的一处热源会逐渐使整个房间变暖。

随机游走与扩散

扩散的微观本质是粒子（如分子）的随机运动。大量粒子的集体随机运动在宏观上表现为扩散现象。在之前的随机游走课程中，我们学到粒子向左或向右移动的概率各为 1/2。考虑一个位置 i 上的粒子浓度（或温度）T_i^n，在下一时间步，该点的粒子可能来自左边的 i-1、右边的 i+1 或留在原地 i。

通过概率分析，可以得到扩散的离散更新公式：

T_i^{n+1} = T_i^n + (D Δt / (Δx)^2) * (T_{i-1}^n - 2T_i^n + T_{i+1}^n)

其中 D 是扩散系数，控制扩散的快慢。公式中 (T_{i-1} - 2T_i + T_{i+1}) 这一项是二阶中心差分，它是空间二阶导数的离散近似。

从离散到连续：扩散方程（热方程）

取 Δt 和 Δx 趋近于零的极限，上述离散方程就变成了著名的扩散方程，也称为热方程：

∂T/∂t = D * ∂²T/∂x²

这里，∂²T/∂x² 是温度在空间上的二阶导数（曲率）。方程表明，某点温度随时间升高的速率，正比于该点温度分布的曲率。在热点（局部极大值，曲率为负），温度会下降；在冷点（局部极小值，曲率为正），温度会上升，从而导致整体趋于平滑。

扩散方程是一种抛物型偏微分方程，它描述的是耗散和 smoothing 过程。

结合平流与扩散 🌪️

在真实的物理系统（如海洋或大气）中，平流和扩散往往是同时发生的。例如，海洋中的热量既会被洋流携带（平流），也会在静止的水层中缓慢传导（扩散）。

组合模型

我们可以简单地通过算子分裂方法来组合这两个过程。即，在一个时间步 Δt 内，先计算纯平流带来的变化，再在平流结果的基础上计算扩散带来的变化（或者交换顺序）。对应的组合离散更新步骤可以写作：

# 平流步（使用改进的格式，如中心差分）
T_inter = advection_step(T^n, u, Δt, Δx)
# 扩散步
T^{n+1} = diffusion_step(T_inter, D, Δt, Δx)

对应的连续偏微分方程就是平流-扩散方程：

∂T/∂t = -u * ∂T/∂x + D * ∂²T/∂x²

这个方程是气候模型中描述能量或物质输运的基础之一。通过数值求解这个方程，我们可以模拟温度、盐度或污染物浓度等在空间中的演化。

总结与展望 📈

本节课我们一起学习了偏微分方程的两个基本构建模块：平流和扩散。

• 平流方程 ∂T/∂t = -u ∂T/∂x 描述了物理量随流体的整体输运。其数值求解需要注意格式选择（如中心差分）和稳定性条件（CFL条件）。
• 扩散方程 ∂T/∂t = D ∂²T/∂x² 描述了物理量由于随机运动导致的从高浓度向低浓度的传播。其离散形式自然引出二阶差分。
• 将两者结合得到的平流-扩散方程 ∂T/∂t = -u ∂T/∂x + D ∂²T/∂x²，能够更真实地模拟许多物理过程，是气候科学和流体力学中的核心方程。

我们从离散的物理图像（质量守恒、随机游走）出发，推导出离散更新公式，再通过取极限得到连续的偏微分方程，最后又回到离散的数值求解。这种“离散-连续-离散”的循环是计算思维在科学计算中的典型体现。

在未来的课程中，我们将把这些概念扩展到二维或三维空间，并应用于更复杂、更真实的气候系统模型中。理解这些基础过程，是构建和理解复杂地球系统模型的关键第一步。

L24：电阻器、模板和气候模型 🧠

在本节课中，我们将学习如何将离散的物理问题（如电阻网络）与连续的数学概念（如偏微分方程）联系起来。我们将通过电阻网络、数值模板（Stencil）的应用，以及一个简化的二维海洋气候模型来探索这些概念。课程的核心是理解拉普拉斯算子（Laplacian）在离散和连续世界中的表现形式。

---

概述 📋

本节课我们将从具体的物理实例——电阻网络出发，推导出描述平衡状态的方程。接着，我们会看到同样的数学形式如何应用于图像处理或求解偏微分方程的数值模板。最后，我们会将这些概念扩展到一个二维的海洋热量传输模型中，这是现代气候模型的一个基础组成部分。

---

从离散到连续：电阻网络与拉普拉斯方程 ⚡

上一节我们概述了课程目标，本节中我们来看看一个具体的物理例子：由1欧姆电阻组成的网络。

想象你手中有一袋1欧姆的电阻器。根据欧姆定律，电压降 V 等于电流 I 乘以电阻 R：V = I * R。对于一个1欧姆的电阻，若有9伏特的电压降，则对应9安培的电流。

我们可以将这些电阻排列成有趣的图案，例如一个“十字形”。在这个图案中，我们关注五个点的电压：中心点电压 V 和其北、东、南、西四个邻居点的电压（V_N, V_E, V_S, V_W）。根据基尔霍夫电流定律，流入一个节点的净电流为零（除非有外部电流注入）。对于中心节点，流经四个电阻的电流之和必须为零。

每个支路的电流等于该支路的电压差除以电阻（1欧姆）。因此，我们可以得到中心节点的方程：
(V_N - V) + (V_E - V) + (V_S - V) + (V_W - V) = I_out
其中 I_out 是可能的外部输出电流。整理后得到：
(V_N + V_E + V_S + V_W) - 4V = I_out

这个方程是理解后续内容的关键。它表明，中心点的电压与其四个邻居点电压的平均值有关，其差值正比于外部电流。

现在，让我们将这个“十字形”网络扩展成一个更大的网格，例如一个包含60个节点（黑色圆点）的网格。每个内部节点都遵循上述同样的规则。如果我们知道所有边界条件（例如某些点接电池或接地）以及可能的外部注入电流，我们就能为每个内部节点列出一个类似的方程，从而得到一个大型线性方程组。

求解这个方程组（例如使用雅可比迭代法或计算机）可以得到网格中每个节点的电压。有趣的是，这个离散的电阻网络问题，在网格无限细化的极限下，会趋近于一个连续的数学问题——拉普拉斯方程或泊松方程。

为了建立这种联系，我们考虑一个定义在整数网格点上的函数 f(x, y)，例如 f(x, y) = 2x^3 - y^4。我们可以计算该函数在中心点 (0,0) 与其邻居点的值。

以下是计算二阶中心差分（离散拉普拉斯算子）的步骤：

1. 计算一阶差分（相邻点的差值）：
   • 北向：f_N - f
   • 南向：f - f_S
   • 东向：f_E - f
   • 西向：f - f_W
2. 计算二阶差分（一阶差分的差分）：
   • 垂直方向二阶差分：(f_N - f) - (f - f_S) = f_N + f_S - 2f
   • 水平方向二阶差分：(f_E - f) - (f - f_W) = f_E + f_W - 2f
3. 将两个方向的二阶差分相加，得到离散拉普拉斯算子：
   (f_N + f_E + f_S + f_W) - 4f

这与我们电阻网络方程的左端形式完全一致！如果我们考虑网格间距为 h（而不一定是1），则需要用 h^2 进行归一化，使其成为导数的近似：
[ (f_N + f_E + f_S + f_W) - 4f ] / h^2 ≈ ∂²f/∂x² + ∂²f/∂y²

右边的连续算子 ∂²f/∂x² + ∂²f/∂y² 被称为拉普拉斯算子，记作 ∇²f 或 Δf。因此，离散的电阻网络方程 (邻居和) - 4*(自身) = I 对应的连续形式就是泊松方程：∇²f = g，其中 g 是已知的源项（对应于离散方程中的 I）。

总结来说，我们展示了如何从一个具体的、离散的物理问题（电阻网络）出发，通过数学推导，连接到描述连续现象的偏微分方程（泊松方程）。这体现了计算思维中离散与连续之间的深刻联系。

---

数值模板：在代码中应用拉普拉斯算子 🎨

上一节我们介绍了拉普拉斯算子在离散和连续形式下的数学表达，本节中我们来看看如何在计算机程序中实现它，这就要用到模板。

模板是一个小型的、带权重的窗口，我们在数据数组（如图像矩阵、温度场）上滑动它。在每一步，将窗口覆盖的数据与模板的权重系数相乘并求和，结果写入输出数组的对应位置。这广泛应用于图像处理（如模糊、边缘检测）和求解偏微分方程。

以我们之前得到的离散拉普拉斯模板为例，它是一个3x3的窗口，中心权重为4，上下左右四个邻居权重为-1，其余角点权重为0。

[ 0, -1,  0 ]
[-1,  4, -1 ]
[ 0, -1,  0 ]

在Julia中实现模板运算时，需要处理边界问题。以下是几种常见的边界处理策略：

• 忽略边界：模板窗口可以滑出数据边界，通常只计算内部点。
• 固定边界：假设边界外的值为某个常数（如0）。
• 周期边界：假设数据在边界处是周期重复的，上边界之外取用下边界的数据，左边界之外取用右边界的数据。
• 零通量边界：在物理模拟中常用，意味着边界没有物质或热量通过。这可以通过设置“幽灵细胞”来实现，使边界上的梯度为零。

在代码中，我们可以利用Julia的两个特性来优雅地实现模板运算：

1. 笛卡尔索引：CartesianIndex 允许我们用一个变量来表示多维索引（如 CartesianIndex(2,3)），方便在循环中统一处理。
2. 偏移数组：OffsetArray 允许数组的索引不从1开始，例如可以从0开始，这便于实现以当前点为中心的模板窗口。

核心代码逻辑是遍历每个内部网格点，收集其邻居窗口的值，与模板权重进行点乘求和：

# 假设 `data` 是输入数组，`stencil` 是模板权重矩阵，`neighborhood` 是相对索引集合
for I in interior_indices
    window = [data[I + offset] for offset in neighborhood]
    result[I] = sum(window .* stencil)
end

通过使用偏移数组预先填充边界外的“幽灵细胞”，我们可以使上述循环代码简洁且无需额外的边界判断，因为所有必要的索引访问都是合法的。

---

迈向气候模型：二维平流-扩散方程 🌊

上一节我们学习了在网格上应用模板的通用技术，本节中我们来看看如何将这些技术应用于一个更复杂、更接近现实的问题：模拟海洋中的热量传输，这是气候模型的核心组件之一。

海洋热量传输主要由两个物理过程驱动：

1. 平流：海水流动（如洋流）将热量从一个地方携带到另一个地方。
2. 扩散：热量从高温区域向低温区域自然扩散。

描述这个过程的连续方程是二维平流-扩散方程：
∂T/∂t = - (u * ∂T/∂x + v * ∂T/∂y) + κ * (∂²T/∂x² + ∂²T/∂y²)
其中：

• T(x, y, t) 是温度场。
• u(x, y) 和 v(x, y) 是流速场的x和y方向分量。
• κ 是热扩散系数。
• - (u * ∂T/∂x + v * ∂T/∂y) 是平流项。
• κ * ∇²T 是扩散项（拉普拉斯算子）。

为了在计算机上求解，我们需要将方程离散化。这涉及到用有限差分来近似所有的偏导数：

• 一阶导数（平流项）可以用中心差分近似：∂T/∂x ≈ (T[i+1, j] - T[i-1, j]) / (2Δx)
• 二阶导数（扩散项）就是我们熟悉的离散拉普拉斯模板：∇²T ≈ (T[i+1,j] + T[i-1,j] + T[i,j+1] + T[i,j-1] - 4T[i,j]) / Δx² （假设 Δx = Δy）

对于边界条件，在海洋模型中，海岸线通常被视为“无通量”边界，即没有热量或海水穿过。这可以通过设置“幽灵细胞”的值，使得边界处的法向梯度为零来实现。例如，在左边界，我们可以设置幽灵细胞的值等于其右侧第一个内部细胞的值。

在编程实现时，良好的软件工程实践是使用抽象。我们可以定义不同的类型来封装不同层次的信息：

struct Grid
    nx::Int
    ny::Int
    Δx::Float64
    Δy::Float64
    x::Vector{Float64}
    y::Vector{Float64}
end

struct OceanModel
    grid::Grid
    u::Matrix{Float64}  # x方向流速场
    v::Matrix{Float64}  # y方向流速场
    κ::Float64
end

这样，主模拟循环可以清晰且模块化。我们可以轻松更换不同的流速场（例如单涡旋场或双涡旋场），而无需重写核心计算逻辑。通过时间步进迭代更新温度场 T，我们就能模拟出热量在海洋中随洋流平流和扩散的动态过程，如上图所示。

---

总结 🎯

本节课中我们一起学习了：

1. 离散与连续的联系：通过分析电阻网络，我们推导出离散拉普拉斯算子，并展示了它在网格细化时如何逼近连续的泊松方程 ∇²f = g。
2. 数值模板的实现：我们探讨了如何在代码中使用模板（Stencil）来应用拉普拉斯等算子，并讨论了处理计算边界的多种策略。
3. 气候模型基础：我们将概念扩展到二维平流-扩散方程，学习了如何离散化并模拟海洋热量传输这一气候科学中的关键过程，并强调了通过抽象（定义Grid、OceanModel等类型）来构建清晰、可维护的科学代码的重要性。

这些内容构成了从具体物理问题到抽象数学方程，再到实际计算机模拟的完整链条，是计算思维在科学计算领域的典型体现。

L25：气候变化建模 🌍

在本节课中，我们将学习如何使用简化的气候模型来理解和应对气候变化问题。我们将重点介绍一个名为MARGOT的模型，它整合了减缓、移除、地球工程和适应四种技术，并通过优化算法来寻找平衡气候目标与经济成本的最佳策略。

---

模型概述与动机

气候变化是我们这个时代面临的重大挑战。本节课介绍的MARGOT模型旨在将复杂的气候问题简化为一个核心、可解释且开源的模型。该模型不仅包含物理气候过程，还整合了经济学和行为响应因素。

MARGOT代表减缓、适应、移除和地球工程优化。其核心思想是将所有关键要素浓缩到一个小型、可解释的模型中，并用Julia语言编写，以便于交互和探索。

如果你在观看视频，可以访问课程网站 computationalthinking.mit.edu，在那里你可以看到我屏幕上的所有可视化内容，并亲自与气候模型进行交互，甚至修改代码。

---

MARGOT模型的核心组件

上一节我们介绍了MARGOT模型的整体目标，本节中我们来看看构成这个模型的四个核心技术。

MARGOT模型主要包含四种应对气候变化的技术：

1. 减缓：减少温室气体排放。在模型中，这通过一个控制变量 m（范围从0到1）来实现，从总排放量 Q 中减去 m * Q。若 m = 1，则表示完全停止排放。

   mitigated_emissions = Q - m * Q
   

2. 移除：从大气中吸收二氧化碳。这可以通过植树等自然方式或直接空气捕获等工程方式实现。该操作会减少大气中的碳存量。

3. 地球工程：改变地球的反射率。例如，向平流层注入微小颗粒以反射阳光，从而为地球降温。这是一种不完美的抵消方法。

4. 适应：调整社会系统以应对已发生的气候变化影响，例如建造海堤或安装空调。这需要付出成本，但可以减少损失。

所有这些控制措施都会产生成本。模型的目标是在减少气候变化损害和避免过度控制成本之间找到最佳平衡。

---

运行模型与初步分析

了解了模型的基本构成后，我们现在可以运行它，并观察不同策略如何影响结果。

在基准情景下，如果不采取任何气候政策，二氧化碳排放将持续增加，导致全球温度上升。例如，将模型中的减缓控制滑块向右移动，意味着我们减少排放，这将直接导致未来温度上升幅度降低。

然而，完全依靠减缓来实现将温升控制在2°C以内的目标（如《巴黎协定》所要求的），可能需要非常激进的措施。这时，我们可以引入第二种技术：碳移除。

通过同时调整减缓和移除，我们开始面临一个有趣的优化问题：哪种组合能以最低成本实现气候目标？因为减缓和移除都有其各自的成本曲线。

---

使用JuMP进行优化

前面我们看到，手动调整策略来寻找最佳组合是困难的。因此，我们需要借助数学优化。在Julia中，我们使用一个名为JuMP的领域特定语言来表述和解决优化问题。

JuMP允许我们以接近数学公式的方式描述问题。例如，假设我们有一个需要最小化的成本函数 cost(x)，并带有变量 x 的约束。

using JuMP, Ipopt

model = Model(Ipopt.Optimizer) # 创建模型，使用Ipopt求解器
@variable(model, -10 <= x <= 10) # 定义变量x及其范围
@NLobjective(model, Min, x^2 + 2) # 定义需要最小化的非线性目标函数
optimize!(model) # 运行优化
value(x) # 获取最优解下x的值

在这个简单例子中，我们寻找函数 f(x) = x^2 + 2 在 x 属于 [-10, 10] 区间内的最小值。优化器会快速找到 x = 0 时取得最小值 2。

对于MARGOT模型，优化问题则复杂得多。我们需要优化的变量是随时间变化的减缓和移除策略曲线，维度高达40维（例如，每10年一个值，共200年）。约束条件包括温度不能超过2°C，以及控制变量每年的变化幅度有上限等。JuMP能够高效地处理这类大规模、带约束的非线性优化问题。

---

探索最优策略与敏感性分析

通过JuMP进行优化后，MARGOT模型可以为我们计算出满足特定温升目标（如2°C）的最优政策组合。结果通常显示，最优路径包含前期积极的减缓和后期逐步增加的碳移除。

模型的一个强大之处在于可以进行敏感性分析。我们可以改变模型的关键假设，观察最优策略如何随之变化：

• 气候敏感性：改变“气候敏感性”参数（即CO2浓度翻倍导致的温升），会影响所需的政策强度。
• 经济贴现率：这个参数反映了我们对未来收益与当前成本的重视程度。贴现率越高，意味着越不关心未来，模型可能更倾向于依赖“创可贴”式的地球工程或未来移除技术，而非当下的深度减排。
• 技术成本假设：如果未来减缓技术成本预期会下降，最优策略可能会包含更多的减缓措施。

这种分析有助于我们理解，最优气候政策并非固定不变，它强烈依赖于我们对科学和经济的底层假设与价值判断。

---

总结与行动建议

本节课中，我们一起学习了如何使用简化的集成评估模型MARGOT来理解和规划气候应对策略。我们看到了减缓、移除、地球工程和适应四种技术的作用，并了解了如何利用JuMP优化包在复杂的约束条件下（如2°C温升目标）寻找成本效益最高的政策路径。

最后，对于有兴趣投身气候领域的同学，建议可以从参与开源气候项目（如MARGOT本身）或加入专注于气候解决方案的初创公司开始，在能源系统、可持续农业等多个方向寻找机会。

气候变化是每个人的问题，而计算思维和开源工具为我们提供了参与解决这一全球性挑战的途径。