Part I · 第 2 章微积分 · 2.5Part I · Ch. 2 Calculus · 2.5

常微分方程Ordinary Differential Equations

微分方程不告诉你"状态是多少"，而是告诉你在每个状态下它会怎样变化。电机转速、关节角度、电池电量、机器人位置，凡是随时间演化的东西都由它描述。解微分方程就像顺着风向图走：从起点出发，每一刻都朝箭头指的方向走一小步。 A differential equation does not tell you "what the state is" but how it changes in each state. Motor speed, joint angle, battery charge, robot position: everything that evolves in time is described by one. Solving it is like following a wind map: start somewhere and at every instant take a small step in the direction the arrow points.

约 30 分钟~30 min 前置：2.1 导数 · 2.3 积分 · 1.11 特征值（推荐）Prereqs: 2.1 · 2.3 · 1.11 Eigenvalues (recommended) 5 个交互5 interactives

1为什么机器人需要它Why robots need this

机器人是一个动力系统：给它一个输入（电压、力矩），它的状态按物理规律随时间演化。写成公式就是 $\dot{\mathbf x}=\mathbf f(\mathbf x,\mathbf u)$：状态的变化率由当前状态和输入决定。

电机的转速、关节的角度、足式机器人的质心位置，都服从微分方程。
控制器的设计，就是挑选输入 $\mathbf u$，让这个方程的解沿着我们想要的方向走，而且稳定。
仿真器（MuJoCo、Isaac、PyBullet）每一步都在数值求解一个微分方程；强化学习里用到的"环境转移"，在连续的世界里就是它的离散化。

这一节的目标：读懂 $\dot x=f(x)$ 的几何意义，会判断一个线性系统稳不稳，理解为什么步长和方法的选择会让仿真失效。

A robot is a dynamical system: give it an input (voltage, torque) and its state evolves in time under physical laws. In symbols $\dot{\mathbf x}=\mathbf f(\mathbf x,\mathbf u)$: the rate of change of the state is fixed by the current state and input.

Motor speed, joint angle and a legged robot's center-of-mass position all obey differential equations.
Controller design means choosing the input $\mathbf u$ so that the solution of this equation goes where we want and stays stable.
Simulators (MuJoCo, Isaac, PyBullet) solve a differential equation numerically at every step; the "environment transition" used in reinforcement learning is, in a continuous world, its discretization.

This lesson's goals: read the geometry of $\dot x=f(x)$, tell whether a linear system is stable, and see why the choice of step size and method can break a simulation.

2直觉：顺着箭头走Intuition: follow the arrows

$\dot x=f(t,x)$ 的意思是：在平面上的每一点 $(t,x)$，都规定了一个斜率。把这些斜率画成一片小线段，就是斜率场（方向场）。解 $x(t)$ 就是一条处处与小线段相切的曲线：它像一片叶子顺着水流漂。换一个起点，就是另一条解，而同一点只会有一条解（在合理的条件下），所以解的曲线不会相交。

在下面的图里，拖动黄点（初始条件），点击空白处再放几条解。换几个方程看看。

$\dot x=f(t,x)$ says: at every point $(t,x)$ of the plane a slope is prescribed. Draw those slopes as little segments and you get a slope field (direction field). A solution $x(t)$ is a curve tangent to the segments everywhere, like a leaf drifting on a current. A different starting point gives a different solution, and through each point there is only one solution (under reasonable conditions), so solution curves never cross.

In the figure below, drag the amber dot (the initial condition) and click empty space to drop more solutions. Try the different equations.

FIG 1斜率场与解Slope field and solutions

小线段是各点的斜率，琥珀色曲线是可拖动初始条件的主解（向前向后都算），青色是点击放下的其他解。注意：曲线一旦"吸附"到某个水平的位置，说明那里是平衡点（斜率为 0）；稳定的平衡点把附近的解都吸引过来。 The small segments are the slopes; the amber curve is the main solution with a draggable initial condition (computed forward and backward), and the cyan ones are the other solutions you click in. When curves settle onto a horizontal level, that level is an equilibrium (slope 0); a stable equilibrium attracts nearby solutions.

选"逻辑斯蒂"，在 $x>0$ 处放下至少 3 条不同起点的解（3 的上方和下方都试试）：它们都靠向同一条水平线（稳定平衡点 $x=3$）。Pick "logistic" and drop at least 3 solutions from different starts with $x>0$ (try both above and below 3): they all approach the same horizontal line (the stable equilibrium $x=3$).
选"$\dot x=0.25x^2$"，把起点放到 $x_0>1$：解在有限时间内爆炸，曲线冲出画面。Pick "$\dot x=0.25x^2$" and put the start at $x_0>1$: the solution blows up in finite time and shoots off the plot.

核心洞见Key insight

微分方程把"随时间怎样变"写成"当前状态的函数"。解不是公式，而是一条被箭头牵着走的轨迹。有了这个图景，"稳定"就是"所有附近的轨迹都靠向这里"，"数值积分"就是"沿着箭头一小步一小步走"。

A differential equation writes "how it changes in time" as "a function of the current state". A solution is not a formula but a trajectory led along by the arrows. With this picture "stable" means "all nearby trajectories flow toward here", and "numerical integration" means "walk along the arrows in small steps".

3动手玩：一阶、二阶、数值方法Play: first order, second order, numerical methods

一阶系统：电机转速First-order system: motor speed

一个带摩擦的电机，转速 $\omega$ 满足 $\tau\,\dot\omega=-\omega+K u$：$\tau$ 是时间常数，$K$ 是稳态增益。输入阶跃 $u=1$，解是 $\omega(t)=K\,(1-e^{-t/\tau})$：每过一个 $\tau$，离终值的距离缩小到原来的 $1/e\approx37\%$；过 $\tau$ 时到达终值的 63%，过 $4\tau$ 时到 98%。输入是正弦时，系统像一个低通滤波器：高频被削弱，还会滞后。

A motor with friction: speed $\omega$ obeys $\tau\,\dot\omega=-\omega+K u$, where $\tau$ is the time constant and $K$ the steady-state gain. For a step input $u=1$ the solution is $\omega(t)=K\,(1-e^{-t/\tau})$: every $\tau$ the distance to the final value shrinks to $1/e\approx37\%$ of what it was; at $t=\tau$ it has reached 63% of the final value, at $4\tau$ 98%. With a sinusoidal input the system acts like a low-pass filter: high frequencies are attenuated and lag behind.

FIG 2一阶响应First-order response

琥珀色是输入 $u(t)$，青色是转速 $\omega(t)$（用 RK4 数值求解）。竖虚线标出 $t=\tau$，横虚线标出终值的 63%。 Amber is the input $u(t)$ and cyan the speed $\omega(t)$ (solved numerically with RK4). The vertical dashed line marks $t=\tau$ and the horizontal one 63% of the final value.

阶跃输入下，调时间常数 $\tau=0.5$ s：转速在 $t=0.5$ s 时到达终值的 63%。With a step input, set $\tau=0.5$ s: the speed reaches 63% of its final value at $t=0.5$ s.
切到"正弦"，调频率使 $\omega_{in}\tau\approx1$：输出幅值降到约 71%，滞后 45°。Switch to "sine" and tune the frequency so $\omega_{in}\tau\approx1$: the output amplitude drops to about 71% and lags by 45°.

二阶系统：质量–弹簧–阻尼器，与极点Second-order system: mass–spring–damper, and poles

关节、悬架、PD 控制的机械臂，常见形式是 $\ddot\theta+2\zeta\omega_n\dot\theta+\omega_n^2\theta=\omega_n^2 r$。$\omega_n$ 是固有频率，$\zeta$ 是阻尼比。把它写成 $\dot{\mathbf x}=A\mathbf x$（状态 $\mathbf x=(\theta,\dot\theta)$），$A$ 的特征值（又叫极点）决定一切：解里含 $e^{\lambda t}$，实部 $\mathrm{Re}\,\lambda<0$ 才会衰减，虚部决定振荡。右图是复平面，拖动紫色的极点看看响应怎么变；极点越靠左越快，越靠近虚轴越慢，跨过虚轴就发散。

Joints, suspensions and PD-controlled arms often take the form $\ddot\theta+2\zeta\omega_n\dot\theta+\omega_n^2\theta=\omega_n^2 r$. $\omega_n$ is the natural frequency and $\zeta$ the damping ratio. Write it as $\dot{\mathbf x}=A\mathbf x$ (state $\mathbf x=(\theta,\dot\theta)$) and the eigenvalues of $A$ (also called poles) decide everything: the solution contains $e^{\lambda t}$, which decays only if $\mathrm{Re}\,\lambda<0$, and the imaginary part sets the oscillation. The right panel is the complex plane; drag the purple pole and watch the response. The further left, the faster; the closer to the imaginary axis, the slower; across it, the response diverges.

FIG 3阻尼比与极点Damping ratio and poles

阶跃响应 θ(t)step response θ(t)

极点（复平面）poles (complex plane)

左：对单位阶跃输入的响应。右：极点位置，紫色点可拖动（$\zeta<1$ 时是一对共轭极点；$\zeta\ge1$ 时两个极点落在实轴上）。灰色半圆是 $|\lambda|=\omega_n$。 Left: the response to a unit step. Right: the pole locations; the purple dot is draggable (a conjugate pair for $\zeta<1$, two real poles for $\zeta\ge1$). The gray arc is $|\lambda|=\omega_n$.

找到临界阻尼 $\zeta=1$：最快且不超调的响应。Find critical damping $\zeta=1$: the fastest response with no overshoot.
让超调量约为 16%（$\zeta\approx0.5$，工程里常用的折中）。Make the overshoot about 16% ($\zeta\approx0.5$, a common engineering compromise).
把阻尼调成负数（拖动极点到虚轴右侧）：响应越摆越大，系统不稳定。Make the damping negative (drag the pole right of the imaginary axis): the response swings ever wider, unstable.

数值方法：步长、方法与"仿真爆炸"Numerical methods: step size, method, and "simulation blow-up"

一个理想的无阻尼振子 $\ddot x=-x$，其能量 $E=\tfrac12(x^2+v^2)$ 应该永远不变，相图里的轨迹是一个圆。用不同的数值方法、不同的步长 $h$ 积分，看谁能保住这个圆。

An ideal undamped oscillator $\ddot x=-x$ has energy $E=\tfrac12(x^2+v^2)$ that should stay constant forever; in the phase portrait the trajectory is a circle. Integrate it with different numerical methods and step sizes $h$ and see which one keeps the circle.

FIG 4谁能保住能量？Who conserves energy?

相图 (x, v)phase portrait (x, v)

log₁₀ (E / E₀) 随时间log₁₀ (E / E₀) versus time

灰色圆是精确解。显式欧拉每步把能量乘上 $(1+h^2)$，轨迹向外螺旋；半隐式（辛）欧拉能量在一个小范围内振荡，永不漂走；RK4 精度高，但在大步长下能量也会缓慢衰减。 The gray circle is the exact solution. Explicit Euler multiplies the energy by $(1+h^2)$ each step and spirals outward; semi-implicit (symplectic) Euler keeps the energy oscillating in a small band and never drifts; RK4 is very accurate, but at large steps its energy also decays slowly.

选显式欧拉，步长 $h\ge0.2$：20 秒内能量涨到初始的 2 倍以上。Pick explicit Euler with $h\ge0.2$: the energy more than doubles within 20 s.
选半隐式欧拉，同样的步长：能量保持在 ±25% 以内，没有漂走。Pick semi-implicit Euler with the same step: the energy stays within ±25% and does not drift.

4再看数学The math

提示：悬停在公式带底色的符号上，FIG 3 中对应的部分会高亮。Tip: hover the highlighted symbols; the matching part of FIG 3 lights up.

标准形式与状态空间Standard form and state space

任何高阶方程都可以写成一阶方程组。例如 $\ddot\theta=g(\theta,\dot\theta,u)$：令 $\mathbf x=(\theta,\dot\theta)$，则 $\dot{\mathbf x}=(\dot\theta,\ g(\theta,\dot\theta,u))$。状态 $\mathbf x$ 是"知道它就能预测未来"的最少信息。

Any higher-order equation can be written as a first-order system. For $\ddot\theta=g(\theta,\dot\theta,u)$, set $\mathbf x=(\theta,\dot\theta)$ to get $\dot{\mathbf x}=(\dot\theta,\ g(\theta,\dot\theta,u))$. The state $\mathbf x$ is the minimum information that lets you predict the future.

$$ \dot{\mathbf x}=\mathbf f(\mathbf x,\mathbf u,t),\qquad \mathbf x(t_0)=\mathbf x_0. $$

只要 $\mathbf f$ 光滑（局部 Lipschitz），给定初值，解存在并且唯一：所以解曲线不相交。但解可能在有限时间内爆炸（FIG 1 的 $\dot x=0.25x^2$）。

As long as $\mathbf f$ is smooth (locally Lipschitz), for a given initial value a solution exists and is unique: hence solution curves never cross. A solution may, however, blow up in finite time (the $\dot x=0.25x^2$ of FIG 1).

线性系统与极点Linear systems and poles

$$ \dot{\mathbf x}=A\mathbf x\ \Longrightarrow\ \mathbf x(t)=e^{At}\mathbf x_0=\sum_i c_i\,e^{\hl{lam}{\lambda_i}t}\,\mathbf v_i,\qquad A\mathbf v_i=\lambda_i\mathbf v_i. $$

沿每个特征向量 $\mathbf v_i$，系统退化成一个一阶方程 $\dot c=\lambda c$，解是 $e^{\lambda t}$（特征值的概念见 1.11 节）。所以：$\mathrm{Re}\,\lambda_i<0$ 对所有 $i$ 成立，则系统渐近稳定；只要有一个 $\mathrm{Re}\,\lambda_i>0$，就不稳定；虚部给出振荡频率。二阶系统 $\ddot\theta+2\zeta\omega_n\dot\theta+\omega_n^2\theta=0$ 的极点是

Along each eigenvector $\mathbf v_i$ the system reduces to a first-order equation $\dot c=\lambda c$ with solution $e^{\lambda t}$ (eigenvalues: Lesson 1.11). So the system is asymptotically stable iff $\mathrm{Re}\,\lambda_i<0$ for all $i$, unstable if some $\mathrm{Re}\,\lambda_i>0$, and the imaginary part gives the oscillation frequency. The poles of $\ddot\theta+2\zeta\omega_n\dot\theta+\omega_n^2\theta=0$ are

$$ \lambda_{1,2}=\hl{sigma}{-\zeta\omega_n}\;\pm\;\hl{wd}{\,j\,\omega_n\sqrt{1-\zeta^2}}\qquad(\zeta<1),\qquad M_p=e^{-\pi\zeta/\sqrt{1-\zeta^2}}. $$

其中 $M_p$ 是阶跃响应的超调量（峰值超出终值的比例）：$\zeta=0.5$ 时约 16%，$\zeta=0.7$ 时约 5%，$\zeta\ge1$ 时为 0。

Here $M_p$ is the step response's overshoot (peak above the final value, as a fraction): about 16% at $\zeta=0.5$, about 5% at $\zeta=0.7$, and 0 for $\zeta\ge1$.

一阶系统 $\tau\dot\omega=-\omega+Ku$ 的极点是 $-1/\tau$：时间常数越小，极点越靠左，响应越快。对非线性系统，在平衡点处用泰勒展开（2.4 节）得到线性化 $\dot{\delta\mathbf x}=J\,\delta\mathbf x$，$J$ 是雅可比，它的特征值给出这个平衡点的局部稳定性。

For the first-order system $\tau\dot\omega=-\omega+Ku$ the pole is $-1/\tau$: the smaller the time constant, the further left the pole and the faster the response. For a nonlinear system, Taylor-expand at an equilibrium (Lesson 2.4) to get the linearization $\dot{\delta\mathbf x}=J\,\delta\mathbf x$; the eigenvalues of the Jacobian $J$ give the equilibrium's local stability.

数值求解Numerical solution

方法Method	一步One step	阶Order
显式欧拉Explicit Euler	$\mathbf x_{k+1}=\mathbf x_k+h\,\mathbf f(\mathbf x_k)$	1
中点法（RK2）Midpoint (RK2)	$\mathbf x_{k+1}=\mathbf x_k+h\,\mathbf f\big(\mathbf x_k+\tfrac h2\mathbf f(\mathbf x_k)\big)$	2
RK4	$\mathbf x_{k+1}=\mathbf x_k+\tfrac h6(\mathbf k_1+2\mathbf k_2+2\mathbf k_3+\mathbf k_4)$	4
半隐式欧拉（辛）Semi-implicit (symplectic) Euler	$v_{k+1}=v_k+h\,a(x_k),\ x_{k+1}=x_k+h\,v_{k+1}$	1

它们都是 2.3 节积分规则的推广（欧拉 = 左矩形，中点法 = 中点规则，RK4 ≈ Simpson）。还有一个 2.3 节没有的问题：稳定性。对 $\dot x=\lambda x$（$\lambda<0$），显式欧拉要求 $h<2/|\lambda|$，否则数值解自己会发散，哪怕真解在衰减。系统里同时存在很快和很慢的模态（"刚性"）时，只能用极小的步长，或改用隐式方法。

They generalize the quadrature rules of Lesson 2.3 (Euler = left rectangle, midpoint = midpoint rule, RK4 ≈ Simpson). One issue Lesson 2.3 did not have is stability. For $\dot x=\lambda x$ ($\lambda<0$) explicit Euler requires $h<2/|\lambda|$; otherwise the numerical solution itself diverges even though the true one decays. When a system has both very fast and very slow modes ("stiffness") you must either use tiny steps or switch to implicit methods.

5机器人中的它：PD 控制的关节In robotics: a PD-controlled joint

一根受重力作用的关节：$J\ddot\theta+b\dot\theta+mg\ell\sin\theta=\tau$。控制器输出力矩 $\tau=k_p(\theta^\star-\theta)-k_d\dot\theta$（PD 控制）。线性化后（并补偿重力时）就是 FIG 3 的二阶系统：$\omega_n=\sqrt{k_p/J}$，$\zeta=\dfrac{k_d+b}{2\sqrt{k_pJ}}$。$k_p$ 像弹簧，决定"多快、多硬"；$k_d$ 像阻尼，抑制振荡。下面让它去追一个 60° 的目标，看重力怎样让它停在错误的位置。

A joint under gravity: $J\ddot\theta+b\dot\theta+mg\ell\sin\theta=\tau$. The controller outputs torque $\tau=k_p(\theta^\star-\theta)-k_d\dot\theta$ (PD control). Linearized (with gravity compensated), this is FIG 3's second-order system: $\omega_n=\sqrt{k_p/J}$, $\zeta=\dfrac{k_d+b}{2\sqrt{k_pJ}}$. $k_p$ acts like a spring, setting "how fast, how stiff"; $k_d$ acts like a damper, suppressing oscillation. Below the joint chases a 60° target; watch how gravity makes it stop at the wrong place.

FIG 5PD 控制的关节A PD-controlled joint

参数：$J=0.5$ kg·m²，$b=0.1$ N·m·s/rad，$mg\ell=4$ N·m。曲线是完整非线性方程的数值解（RK4），灰色带是 ±2% 的目标范围。 Parameters: $J=0.5$ kg·m², $b=0.1$ N·m·s/rad, $mg\ell=4$ N·m. The curve is the numerical solution of the full nonlinear equation (RK4); the gray band is ±2% around the target.

不补偿重力，$k_p\le3$：关节停在离目标超过 10° 的地方（稳态误差）。再打开重力补偿，误差消失。Without gravity compensation and $k_p\le3$: the joint stops more than 10° short (steady-state error). Turn gravity compensation on and the error disappears.
打开重力补偿，调 $k_p$、$k_d$，使超调 小于 5%，且 1.5 秒内进入 ±2% 范围。With gravity compensation on, tune $k_p$ and $k_d$ for overshoot below 5% and settling inside ±2% within 1.5 s.

真实系统里In real systems

真实控制器是在离散时间里运行的：每个周期（如 1 kHz）读传感器、算力矩、写给电机。这相当于对上面的连续方程做了一次"零阶保持"的离散化；$k_p$ 太大、周期太长，闭环极点会越过离散系统的稳定边界，机械臂会开始高频抖动。所以控制器的采样周期要比时间常数小得多，这和上面欧拉法的步长限制是同一件事。后续 Part IV 会系统地讲 PID、LQR 和 MPC。

A real controller runs in discrete time: every cycle (say 1 kHz) it reads sensors, computes torque and writes to the motor. That is a "zero-order hold" discretization of the continuous equation above; if $k_p$ is too large for the cycle time, the closed-loop poles cross the discrete stability boundary and the arm starts to chatter at high frequency. So the control period must be much shorter than the time constants, the same restriction as Euler's step limit above. Part IV later covers PID, LQR and MPC systematically.

6常见误区Common pitfalls

误区 1：步长越大越快，所以用大步长Pitfall 1: bigger steps are faster, so use big steps

步长必须小于系统最快模态的时间尺度。显式欧拉在 $h>2/|\lambda|$ 时发散，而且常常不是"算不准"，是"炸掉"：数值变成 NaN。仿真里突然出现的"机器人飞走了"，很多就是这个。

The step must be below the time scale of the system's fastest mode. Explicit Euler diverges for $h>2/|\lambda|$, and often not as "inaccurate" but as "blown up": the numbers turn NaN. Many sudden "the robot flew away" moments in simulation are exactly this.

误区 2：欧拉法对振荡系统也够用Pitfall 2: Euler is good enough for oscillators

显式欧拉每一步都往能量里"注入"一点，纯振荡的系统会越摆越大。刚体动力学里常用半隐式欧拉或更高阶方法。

Explicit Euler injects a little energy each step, so a pure oscillator swings ever wider. Rigid-body dynamics typically uses semi-implicit Euler or higher-order methods.

误区 3：只看稳不稳，不看快不快、抖不抖Pitfall 3: looking only at stability, not speed or ringing

所有极点都在左半平面只保证最终会收敛。靠近虚轴的极点收敛得极慢；虚部很大的极点收敛时伴随强烈的振荡。设计时要同时看"实部（快慢）"和"阻尼比（抖不抖）"。

All poles in the left half-plane only guarantee eventual convergence. A pole near the imaginary axis converges very slowly; a pole with a large imaginary part rings strongly. Design looks at both the real part (speed) and the damping ratio (ringing).

误区 4：线性化的稳定性适用于一切Pitfall 4: linearized stability applies everywhere

线性化只在平衡点附近有效（2.4 节）。FIG 5 里不补偿重力时，系统停在的并不是目标，而是另一个平衡点。非线性系统可以有多个平衡点，每个都有自己的局部稳定性。

Linearization is valid only near the equilibrium (Lesson 2.4). In FIG 5 without gravity compensation the joint settles not at the target but at a different equilibrium. A nonlinear system can have several equilibria, each with its own local stability.

7小练习Check yourself

1. $\dot x=-2x$，$x(0)=3$。$x(1)$ 约是多少？1. $\dot x=-2x$, $x(0)=3$. About what is $x(1)$?

解是 $x(t)=3e^{-2t}$，$x(1)=3e^{-2}\approx0.406$。The solution is $x(t)=3e^{-2t}$, so $x(1)=3e^{-2}\approx0.406$.

2. 一阶系统时间常数 $\tau=0.2$ s，对单位阶跃输入，$t=0.2$ s 时输出大约到终值的多少？2. A first-order system with $\tau=0.2$ s gets a unit step. At $t=0.2$ s the output has reached about what fraction of its final value?

$1-e^{-1}\approx0.632$。到 $2\tau$ 时是 86%，到 $4\tau$ 时是 98%。$1-e^{-1}\approx0.632$. At $2\tau$ it is 86% and at $4\tau$ 98%.

3. 线性系统的极点是 $-1\pm2j$。它的行为是？3. A linear system has poles $-1\pm2j$. It behaves:

实部 −1 < 0：稳定，衰减的包络是 $e^{-t}$；虚部 ±2：振荡角频率是 2 rad/s。Real part −1 < 0: stable with decay envelope $e^{-t}$; imaginary part ±2: oscillation at 2 rad/s.

4. 二阶系统阻尼比 $\zeta=1$（临界阻尼），阶跃响应是？4. A second-order system with damping ratio $\zeta=1$ (critical damping). The step response is:

两个极点重合在 $-\omega_n$。$\zeta<1$ 会超调，$\zeta>1$ 的慢极点更靠近原点，响应更慢。Both poles coincide at $-\omega_n$. $\zeta<1$ overshoots and for $\zeta>1$ the slow pole is nearer the origin, so the response is slower.

5. 用显式欧拉模拟一个无阻尼振子，时间一长会看到什么？5. You simulate an undamped oscillator with explicit Euler. After a long time you see:

显式欧拉每步把能量乘 $1+h^2$，指数增长。半隐式欧拉则没有这个漂移。Explicit Euler multiplies the energy by $1+h^2$ per step: exponential growth. Semi-implicit Euler has no such drift.

8延伸阅读Further reading

3Blue1Brown, Differential Equations — 第 1 章、第 2 章：相空间与方向场Ch. 1 and 2: phase space and direction fields
S. Strogatz, Nonlinear Dynamics and Chaos — 一维、二维流与稳定性one- and two-dimensional flows and stability
K. Åström & R. Murray, Feedback Systems — 第 3、5 章：动态行为与线性系统Ch. 3 and 5: dynamic behavior and linear systems
E. Hairer, S. Nørsett & G. Wanner, Solving Ordinary Differential Equations I — 数值方法的经典教材the classic reference on numerical methods
K. Lynch & F. Park, Modern Robotics, Ch. 11 — 机器人控制：二阶系统与 PD 控制robot control: second-order systems and PD control

方法Method	一步One step	阶Order
显式欧拉Explicit Euler	\(\mathbf x_{k+1}=\mathbf x_k+h\,\mathbf f(\mathbf x_k)\)	1
中点法（RK2）Midpoint (RK2)	\(\mathbf x_{k+1}=\mathbf x_k+h\,\mathbf f\big(\mathbf x_k+\tfrac h2\mathbf f(\mathbf x_k)\big)\)	2
RK4	\(\mathbf x_{k+1}=\mathbf x_k+\tfrac h6(\mathbf k_1+2\mathbf k_2+2\mathbf k_3+\mathbf k_4)\)	4
半隐式欧拉（辛）Semi-implicit (symplectic) Euler	\(v_{k+1}=v_k+h\,a(x_k),\ x_{k+1}=x_k+h\,v_{k+1}\)	1