🔬
动力学的数值计算方法
欧拉显式 · 欧拉隐式 · 中点差分法 · 蛙跳法
——谁更"高效"？谁更"稳定"？
高中一年级 · 物理计算思维
📖 课程内容
1
预备知识
先看导数、三角函数、振动方程和相空间这些必备概念
2
为什么需要数值方法？
当方程太难算时，计算机怎样一步一步逼近正确答案
4
蛙跳法
位置和速度错开半步更新，长期更稳，更适合振动和轨道问题
5
误差与能量守恒
四种方法到底差多少？用弹簧振子实验来对比
💡
本教程会用到弹簧振子作为贯穿例子。弹簧被压缩后来回振动——这是高中物理里最经典的运动之一！
第 1 / 10 页
📍 先认识导数记号：一个点、两个点是什么意思？
在物理里，我们常用位置 $x(t)$ 表示物体在时刻 $t$ 的位置。只要位置会随时间变化，就可以研究它“变化得有多快”。高中里最自然的想法，是先看一小段时间内的平均变化率。
$$\frac{\Delta x}{\Delta t}=\frac{x(t+\Delta t)-x(t)}{\Delta t}$$
当 $\Delta t$ 越来越小，平均变化率就逼近某一时刻的瞬时变化率。也就是说，当时间取得越来越靠近某个时刻时，变化率也会越来越靠近某个确定的值，这就是导数：
$$\dot{x}=\frac{dx}{dt}=\lim_{\Delta t\to 0}\frac{x(t+\Delta t)-x(t)}{\Delta t}$$
这时，$\dot{x}$ 读作“x 一点”，表示位置对时间的一阶导数，也就是速度。如果再对时间求一次导数，就得到：
$$\ddot{x}=\frac{d^2x}{dt^2}=\frac{d}{dt}\left(\dot{x}\right)$$
它读作“x 两点”，表示速度对时间的变化率，也就是加速度。所以后面方程里的“点”，本质上是在表示“对时间求导”。
这里有三个最关键的对应关系：
$x(t)$：位置
描述物体在每个时刻处在什么位置，是运动问题最基础的量。
$\dot{x}$：速度
表示位置变化得有多快，是由平均变化率取极限得到的瞬时变化率。
$\ddot{x}$：加速度
表示速度变化得有多快，也就是对时间再求一次导数。
💡
所以后面看到 $\dot{x}$ 不要慌，它就是速度；看到 $\ddot{x}$，它就是加速度；看到 $\ddot{x}=-\omega^2x$，意思就是“加速度由当前位置决定”。
这里对“极限”的直观表述，参考了 OpenStax 的教材思路：当自变量越来越接近某个值时，函数值也越来越接近并保持靠近某个确定值。
🌸 从弹簧振子说起
一个弹簧连着小球，被压缩后释放，小球会来回振动。物理上这个运动满足：
$$m\ddot{x} = -kx \quad \Longrightarrow \quad \ddot{x} = -\omega^2 x, \quad \omega = \sqrt{\frac{k}{m}}$$
这个方程有精确解：$x(t) = A\cos(\omega t + \varphi)$。
但现实中的方程往往复杂得多（比如空气阻力、电磁力同时存在），根本无法写出公式。这时候，我们用计算机"一小步一小步地逼近"答案——这就是数值方法。
📐 为什么解析解里会出现余弦？
初中学习三角函数时，角的范围通常只限于锐角。到了高中，三角函数的定义会推广到任意角，也就是角度可以不断增大、不断旋转。
$$\cos\theta=\frac{x(t)}{r}$$
这里 $r$ 是点到原点的距离，$x(t)$ 是转过角 $\theta$ 后，对应点在横轴上的坐标，所以 $x(t)$ 的范围总在 $-r$ 到 $r$ 之间。如果取 $r=1$，就得到单位圆上的余弦函数，此时：
$$-1 \le \cos\theta \le 1$$
当角 $\theta$ 不断变化时，点在圆上转动，横坐标 $x$ 就会周期性地在最大值和最小值之间来回变化。因此，余弦函数天然适合描述往复振动。
$$x(t)=A\cos(\omega t+\varphi)$$
在这个式子里，$\omega t+\varphi$ 叫作三角函数的宗量，它是一个实数，会随着时间不断变化；于是余弦值也会周期性变化。前面的 $A$ 决定最大振幅，$\omega$ 决定变化快慢，$\varphi$ 决定起始相位。
💡
所以弹簧振子的解析解写成余弦，不是“凭空猜的”，而是因为余弦函数本身就具有“周期来回变化”的性质。
第 2 / 10 页
🧠 哈密顿是谁？什么叫哈密顿系统？
威廉·罗恩·哈密顿（William Rowan Hamilton）是 19 世纪爱尔兰数学家、物理学家。他对光学、力学和代数都做出过重要贡献。力学中常说的哈密顿力学，就是以他的名字命名的。
在牛顿力学里，我们常用公式
$$F=ma$$
来描述“力怎样决定加速度”；而在哈密顿力学里，我们更喜欢用位置和动量来描述系统，并关注系统的能量如何决定整个运动：
$$H=H(q,p)$$
这里 $q$ 表示位置，$p$ 表示动量，$H$ 叫作哈密顿量，很多情况下它就可以理解成系统的总能量。如果一个系统的运动规律可以由这样一个能量函数来统一决定，我们就把它叫作哈密顿系统。
看什么量？
哈密顿系统重点看位置 $q$ 和动量 $p$，而不只是看位置和加速度。
由什么决定运动？
由哈密顿量 $H$ 决定。对很多保守系统来说，$H$ 就可以直接理解成总能量。
为什么和这里有关？
弹簧振子正是最简单的哈密顿系统之一，所以比较数值方法时，我们特别关心它们是否尊重总能量和相空间轨道。
🔁 由牛顿方程到哈密顿动力学方程
以一维弹簧振子为例，牛顿方程写成：
$$m\ddot{x}=-kx$$
这是一个二阶微分方程，因为里面出现了二阶导数 $\ddot{x}$。现在引入动量：
$$p=m\dot{x}$$
同时把哈密顿量写成：
$$H(x,p)=\frac{p^2}{2m}+\frac{kx^2}{2}$$
这时，运动规律可以改写成两个一阶方程：
$$\dot{x}=\frac{\partial H}{\partial p}=\frac{p}{m}$$
$$\dot{p}=-\frac{\partial H}{\partial x}=-kx$$
📘
这里出现的“对 $H$ 求偏导数”已经超出了高中阶段的知识范围，可以先把它当成大学数学中的一种运算规则。高中阶段更重要的是看懂最后得到的物理关系。
$$\dot{x}=\frac{p}{m}=v$$
$$\dot{p}=F=-kx$$
这两条在高中阶段是完全可以理解的：第一条表示“位置的一阶导数就是速度”；第二条表示“动量的变化率等于力”，而弹簧振子的回复力正是 $F=-kx$。
✅
所以哈密顿方程中的两个方程，可以这样理解：第一个是速度方程，第二个本质上就是牛顿第二定律在这个系统中的写法。
发生了什么变化？
原来一个二阶方程，被改写成了两个一阶方程。方程的阶数下降了。
代价是什么？
代价是方程的个数增加了一倍：从 1 个方程变成了 2 个方程。
好处是什么？
位置 $x$ 和动量 $p$ 的地位变得对等了。它们一起构成相空间中一个代表系统状态的“坐标点”。
💡
这就是哈密顿力学非常有力量的地方：它把“运动方程”变成了“相空间中一个点如何随时间演化”的问题。
🪐 哈密顿力学视角：我们在追踪“总能量轨道”
哈密顿力学一个非常重要的思想是：把原来“随时间变化的动力学问题”，转化成相空间里的几何问题来研究。
📌
相空间就是用“位置”和“动量”作为坐标轴组成的空间。对一维问题来说，相空间中的一个点就表示系统在某一时刻的完整状态。
对简谐振子来说，它不只是“位置随时间变化”的问题，它还是一个最典型的哈密顿系统。在这个视角下，系统由位置 $x$ 和动量 $p=mv$ 一起决定，总能量写成：
$$H(x,p)=\frac{p^2}{2m}+\frac{kx^2}{2}$$
这里的 $H$ 就是哈密顿量，也就是我们熟悉的总能量。如果系统没有摩擦，真实运动应该一直沿着同一条等能量轨道运动，而不是越绕越大，也不是越绕越小。
$$\frac{p^2}{2m}+\frac{kx^2}{2}=E$$
这就是相空间中的一条椭圆。也就是说：对弹簧振子而言，系统随时间演化，在相空间里看起来就是沿着一条椭圆轨道不断绕行。
真实系统：在相空间里沿椭圆轨道运动
在相空间里，简谐振子的真实轨道是一条椭圆。椭圆的大小对应总能量，轨道不应自己膨胀或收缩。
这正是后面比较数值方法的关键
如果某种方法把椭圆轨道算得越来越外，说明它在“假增能”；如果越来越里，说明它在“假耗能”。更好的方法不只要局部算得准，还要长期尊重这条相空间中的椭圆轨道。
总能量 = 哈密顿量
相空间椭圆
长期结构保真
所以这份课件其实也在回答一个问题：哪种算法更“尊重”哈密顿系统在相空间中的几何结构？
🧭 核心思路：把时间切成小片
把整段时间 $[0, T]$ 切成 $N$ 个小间隔，每段长度叫时间步长 $\Delta t$（也写作 $dt$）。
$$t_0 = 0,\quad t_1 = \Delta t,\quad t_2 = 2\Delta t,\quad \ldots,\quad t_N = T$$
用 $x_n$ 表示在 $t_n$ 时刻的位置，$v_n$ 表示速度。每次从 $(x_n, v_n)$ 推算下一步 $(x_{n+1}, v_{n+1})$，就像下楼梯——一步一步往前走。
🏃 → 🏃 → 🏃 → 🏃 → 🎯
每一步用上一步的信息推算，一步一步逼近真实轨迹
📌 关键公式：泰勒展开的直觉
如果已知现在的状态，下一小步的状态可以近似为：
$$x(t + \Delta t) \approx x(t) + \dot{x}(t)\cdot\Delta t + \underbrace{\frac{\ddot{x}(t)}{2}\Delta t^2 + \cdots}_{\text{误差项}}$$
不同方法的区别，就在于如何取近似、用哪个时刻的信息。
🖼️ 图像直觉：用切线去追曲线
真实运动是一条弯曲轨迹
如果直接知道整条曲线当然最好，但多数复杂方程做不到，所以只能抓住“此刻的变化趋势”去猜下一小步。
数值方法就是连续做“小步预测”
步长越小，每一小段切线预测通常越贴近真实曲线。不同算法的本质区别，就是每一步用哪个方向、方向估得准不准。
第 3 / 10 页
🚀 核心思想：大胆往前冲
用 当前时刻的速度和加速度，直接外推下一步的位置和速度。
📐 递推公式
已知第 $n$ 步的位置 $x_n$ 和速度 $v_n$，计算加速度 $a_n = -\omega^2 x_n$，然后：
$$\boxed{x_{n+1} = x_n + v_n \cdot \Delta t}$$
$$\boxed{v_{n+1} = v_n + a_n \cdot \Delta t = v_n - \omega^2 x_n \cdot \Delta t}$$
$$\boxed{
\begin{bmatrix}x_{n+1}\\ v_{n+1}\end{bmatrix}
=
\begin{bmatrix}
1 & \Delta t\\
-\omega^2\Delta t & 1
\end{bmatrix}
\begin{bmatrix}x_n\\ v_n\end{bmatrix}
}$$
注意：$a_n$ 用的是 $x_n$（当前位置），这就是"显式"的含义——新值可以直接算出来，不需要解方程。
欧拉显式法：用当前点的斜率沿切线方向外推下一小步。
🔢 手算一步（$\omega=1,\ \Delta t=0.1$）
初始条件：$x_0 = 1.0,\ v_0 = 0.0$
1
计算加速度
$a_0 = -\omega^2 x_0 = -1^2 \times 1.0 = -1.0$
2
更新位置
$x_1 = x_0 + v_0 \cdot \Delta t = 1.0 + 0.0 \times 0.1 = 1.0$
3
更新速度
$v_1 = v_0 + a_0 \cdot \Delta t = 0.0 + (-1.0) \times 0.1 = -0.1$
精确解：$x_1^{\text{exact}} = \cos(0.1) \approx 0.9950$，误差 $\approx 0.005$
⚖️ 优缺点
✅ 优点
- 公式极其简单
- 每步计算量最少
- 编程实现只需 2 行
❌ 缺点
- 误差随时间累积
- 能量会持续增大（振幅漂移）
- 步长稍大就发散
⚠️
程序测试中，欧拉显式在 $dt=0.05$、$T=10\text{s}$ 条件下，能量漂移高达 +64.36%！弹簧越振越猛，这显然是不对的。
🧭 图形理解：为什么会“越走越外”
每一步都沿“当前切线”冲出去
对振动问题来说，真实轨迹应当沿着封闭圆轨道转动。但显式欧拉总是用旧信息外推，容易一步一步偏到圆外面去。
为什么偏外就代表能量增大？
在简谐振子的相空间中，离原点越远，代表总能量越大。显式欧拉走出的点半径不断变大，所以会表现成“振幅越来越大、能量越积越多”。
当前位置决定加速度
旧斜率外推
误差持续同向积累
这类方法适合快速试算，但不适合长时间模拟保守系统。
第 4 / 10 页
🛡️ 核心思想：先想后步
用 下一时刻的加速度来更新速度，而不是当前时刻的。但"下一刻的加速度"还不知道——所以需要解个小方程。
📐 递推公式推导
隐式欧拉的定义：
$$x_{n+1} = x_n + v_{n+1} \cdot \Delta t$$
$$v_{n+1} = v_n + a(x_{n+1}) \cdot \Delta t = v_n - \omega^2 x_{n+1} \cdot \Delta t$$
两式联立，把 $v_{n+1}$ 代入第一式：
$$x_{n+1} = x_n + \bigl(v_n - \omega^2 x_{n+1} \Delta t\bigr)\Delta t$$
$$x_{n+1}\bigl(1 + \omega^2 \Delta t^2\bigr) = x_n + v_n \Delta t$$
$$\boxed{x_{n+1} = \frac{x_n + v_n \Delta t}{1 + \omega^2 \Delta t^2}}$$
$$\boxed{v_{n+1} = v_n - \omega^2 x_{n+1} \cdot \Delta t}$$
$$\boxed{
\begin{bmatrix}x_{n+1}\\ v_{n+1}\end{bmatrix}
=
\frac{1}{1+\omega^2\Delta t^2}
\begin{bmatrix}
1 & \Delta t\\
-\omega^2\Delta t & 1
\end{bmatrix}
\begin{bmatrix}x_n\\ v_n\end{bmatrix}
}$$
虽然多了一步代数，但公式仍然可以直接算——这是简谐振子的幸运之处。
欧拉隐式法：用下一时刻的信息反过来校正这一小步。
🔢 手算一步（$\omega=1,\ \Delta t=0.1$）
初始条件：$x_0 = 1.0,\ v_0 = 0.0$
1
计算下一步位置
$x_1 = \dfrac{1.0 + 0.0 \times 0.1}{1 + 1^2 \times 0.1^2} = \dfrac{1.0}{1.01} \approx 0.9901$
2
计算下一步速度
$v_1 = 0.0 - 1^2 \times 0.9901 \times 0.1 = -0.0990$
精确解：$x_1^{\text{exact}} \approx 0.9950$，误差 $\approx 0.005$（量级相同）
⚖️ 优缺点
❌ 缺点
- 能量会持续减小（耗散）
- 振幅越来越小（失真）
- 复杂方程需迭代求解
⚠️
程序测试中，欧拉隐式能量漂移为 −39.16%！弹簧越振越弱，就好像装了减震器——也是失真的。
🧭 图形理解：为什么会“越走越里”
先看“未来位置”，会更保守
隐式欧拉相当于先问一句：“如果我已经到下一步了，那里的加速度是多少？”这种做法会把点往内侧拉，带来更强的稳定性。
稳定的代价是数值耗散
它不容易炸掉，但会把本来应当保存的能量一点点“吃掉”。图上表现为轨道半径不断缩小，物理上就像系统里凭空出现了阻尼。
使用未来信息
数值稳定强
振幅会被压小
因此它很适合“稳定比保真更重要”的刚性问题。
第 5 / 10 页
🎯 核心思想：先试探，再决策
先用欧拉显式走"半步"得到一个预测的中点状态，再用中点的信息推算整步。这样用的是"中间时刻"的斜率，比只用起点更准确！
📐 递推公式（预测-校正）
第一步：预测中点（走半步）
$$x_{\text{mid}} = x_n + v_n \cdot \frac{\Delta t}{2}$$
$$v_{\text{mid}} = v_n + a(x_n) \cdot \frac{\Delta t}{2}$$
第二步：用中点信息走完整步
$$\boxed{x_{n+1} = x_n + v_{\text{mid}} \cdot \Delta t}$$
$$\boxed{v_{n+1} = v_n + a(x_{\text{mid}}) \cdot \Delta t}$$
$$\boxed{
\begin{bmatrix}x_{n+1}\\ v_{n+1}\end{bmatrix}
=
\begin{bmatrix}
1-\frac12\omega^2\Delta t^2 & \Delta t\\
-\omega^2\Delta t & 1-\frac12\omega^2\Delta t^2
\end{bmatrix}
\begin{bmatrix}x_n\\ v_n\end{bmatrix}
}$$
其中 $a(x_{\text{mid}}) = -\omega^2 x_{\text{mid}}$。
中点法：先预测半步，再用中点斜率完成整步推进。
🔢 手算一步（$\omega=1,\ \Delta t=0.1$）
初始条件：$x_0 = 1.0,\ v_0 = 0.0$
1
计算中点位置和速度
$x_{\text{mid}} = 1.0 + 0.0 \times 0.05 = 1.0$
$v_{\text{mid}} = 0.0 + (-1.0) \times 0.05 = -0.05$
2
用中点推算下一步位置
$x_1 = 1.0 + (-0.05) \times 0.1 = 0.9950$
3
用中点推算下一步速度
$a_{\text{mid}} = -1^2 \times 1.0 = -1.0$
$v_1 = 0.0 + (-1.0) \times 0.1 = -0.1$
精确解：$x_1^{\text{exact}} = \cos(0.1) \approx 0.9950$，误差 $\approx 0.000005$。在这个“一步演示”的例子里，它比欧拉法小了大约 1000 倍。
⚖️ 优缺点
✅ 优点
- 精度比欧拉法高一个量级
- 能量近似守恒（仅漂移 0.03%）
- 适合振动、轨道等保守系统
❌ 缺点
- 每步多算一次加速度
- 实现稍复杂
- 极端步长仍可能不稳定
🧭 图形理解：为什么“半步”会更准
真正重要的是“中间方向”
如果曲线在这一小段里明显弯曲，只用起点斜率往前推会偏得比较多。先走半步后再看方向，就能得到更接近整段平均变化趋势的斜率。
这相当于“先试探，再正式前进”
中点法不是盲目一步到位，而是先用一个便宜的预测，摸清这半段路的大致弯曲方向，再用更合理的斜率完成整步推进。
半步预测
中点斜率更代表整体
误差阶数提升
所以它在轨道、振动一类问题里，通常会比欧拉法明显更可信。
欧拉隐式
先假设自己到终点，再用终点信息回头修正。
蛙跳法
让速度落在半步时刻，再和位置交替更新，长期轨道更不易变形。
第 6 / 10 页
🐸 核心思想：位置走整步，速度走半步
蛙跳法不把位置和速度都放在同一个时刻更新，而是让速度落在“半步时刻”。这种交错更新特别适合简谐振子、轨道这类保守系统。
📐 递推公式：交错半步更新
核心递推关系是：
$$v_{n+\frac12} = v_{n-\frac12} + a(x_n)\Delta t$$
$$x_{n+1} = x_n + v_{n+\frac12}\Delta t$$
$$\boxed{
\begin{bmatrix}x_{n+1}\\ v_{n+\frac12}\end{bmatrix}
=
\begin{bmatrix}
1-\omega^2\Delta t^2 & \Delta t\\
-\omega^2\Delta t & 1
\end{bmatrix}
\begin{bmatrix}x_n\\ v_{n-\frac12}\end{bmatrix}
}$$
在代码里，为了启动它，会先用一次半步欧拉得到初始半步速度：
$$v_{\frac12}=v_0+a(x_0)\frac{\Delta t}{2}$$
蛙跳法：位置和速度交错半步，轮流向前推进。
时间轴视角：位置在整数时刻，速度在半整数时刻。
💡
它和中点法一样是二阶方法，但更突出的优点不是“单步最精确”，而是长期保持相空间轨道形状更好。
🧭 图形理解：为什么它长期更稳
不会持续向外发散，也不会持续向内收缩
蛙跳法常见的表现不是单向漂移，而是在真实闭合轨道附近小幅摆动，所以长期看更像真实系统。
为什么这点重要？
对振动和轨道问题，长期结构是否保真，往往比某一步局部误差小一点更重要。这正是蛙跳法最有价值的地方。
二阶方法
长期能量有界振荡
适合保守系统
⚖️ 优缺点
✅ 优点
- 二阶精度，长期表现通常优于欧拉法
- 能量通常不会单向漂移，而是在常值附近小幅振荡
- 适合弹簧振动、轨道、分子动力学等保守系统
❌ 缺点
- 速度和位置不在同一时刻，解释上更绕
- 实现比欧拉法复杂
- 不是所有耗散型问题的首选
第 7 / 10 页
🧪 先把两个模型分清楚
🧭
前面课堂里的四种方法比较，使用的是
简谐振子模型；顶部下载的
Dynamics 演示代码 当前配置演示的是
盒中粒子在重力和阻尼作用下的运动。这两个模型都能用来讲数值积分，但物理问题并不相同，所以这里把它们拆成“课堂模型交互演示”和“真实代码输出回放”两部分来展示。
🎛️ 课堂模型交互演示：浏览器现场算简谐振子
这里不是播放录制好的图片，而是根据你设定的参数，直接在前端按四种递推公式重新计算整段轨迹。
相空间轨迹
横轴是位置 $x$，纵轴是速度 $v$。对简谐振子来说，越偏离闭合轨道，代表数值方法越容易改变系统能量。
📊 截断误差：每步差多少？
每走一步，由于我们舍弃了高阶项，产生的误差叫局部截断误差：
$$\text{欧拉显式/隐式：} \quad \varepsilon_{\text{local}} \sim O(\Delta t^2)$$
$$\text{中点法/蛙跳法：} \quad \varepsilon_{\text{local}} \sim O(\Delta t^3)$$
走了 $N = T/\Delta t$ 步之后，误差累积为全局截断误差：
$$\text{欧拉显式/隐式（一阶方法）：} \quad \varepsilon_{\text{global}} \sim O(\Delta t)$$
$$\text{中点法/蛙跳法（二阶方法）：} \quad \varepsilon_{\text{global}} \sim O(\Delta t^2)$$
💡
步长缩小为原来的 $\frac{1}{2}$：欧拉法误差减半；中点法和蛙跳法这类二阶方法，误差通常减小为 $\frac{1}{4}$。
⚡ 能量漂移（根据当前控制条件实时计算）
下面的百分比会随着你上面的参数一起变化。正值表示总能量偏大，负值表示总能量被数值方法“吃掉”了。
📈 能量随时间变化示意
━━ 欧拉显式（发散）
━━ 欧拉隐式（耗散）
━━ 中点法（近似守恒）
━━ 蛙跳法（有界振荡）
🔍 从物理角度理解
| 方法 | 能量行为 | 物理比喻 |
|---|
| 欧拉显式 | 能量持续增大 | 弹簧越振越猛，像被“加力”了 |
| 欧拉隐式 | 能量持续减小 | 弹簧越振越弱，像装了阻尼 |
| 中点法 | 近似守恒 | 总体接近真实弹簧振动 |
| 蛙跳法 | 有界振荡 | 像在真实轨道附近轻微摆动 |
🎬 Dynamics 真实输出回放：读取 examples/case01/output/display.json
下面这块直接读取 codes/Dynamics/examples/case01/output/display.json。数据来源于 examples/case01/output/display.txt 导出的真实结果，而不是手写示意曲线。
盒中粒子回放
正在加载 Dynamics 输出...
第 8 / 10 页
❓ 第 1 题
欧拉显式方法在模拟弹簧振动时，能量会如何变化？
❓ 第 2 题
中点差分法为什么比欧拉法精度更高？
❓ 第 3 题
如果你在做化工反应仿真，方程变化非常剧烈（刚性方程），应该选哪种方法？
❓ 第 4 题
时间步长 $\Delta t$ 从 $0.1$ 减小到 $0.05$，对于二阶方法（中点法或蛙跳法），全局误差通常变为原来的多少？
🏆
恭喜你完成学习！
你已经掌握了数值积分四种方法的核心原理，这是计算物理的重要基础。
继续探索，物理的世界在等你！
第 10 / 10 页