从游戏物理到金融模型：微分方程数值求解实战指南

1. 项目概述：当物理引擎遇上金融模型

如果你玩过《愤怒的小鸟》或者《荒野大镖客》，一定对里面逼真的物体碰撞、布料飘动或者水波荡漾的效果印象深刻。这些效果背后，是游戏物理引擎在实时求解一系列复杂的微分方程。而当你打开股票软件，看到那些预测股价走势的曲线模型时，背后同样是微分方程在发挥作用。一个在虚拟世界里模拟炮弹轨迹，一个在现实世界中预测资产价格，看似风马牛不相及的两个领域，其数学核心却惊人地一致——微分方程。

这个项目，就是一次打破次元壁的实践。我们将从最直观、最好玩的游戏开发工具Unity入手，亲手搭建一个物理世界，感受微分方程如何驱动小球滚动、弹簧振动。然后，我们会带着这份对微分方程的“手感”，转向数据分析的利器Python，去构建一个简单的金融模型，比如模拟股票价格的随机游走。我的目的不是把你培养成数学博士，而是让你通过“做”来理解微分方程这个抽象概念。你会发现，那些令人生畏的数学符号，在代码的世界里，不过是一行行控制物体行为或数据演变的指令。无论你是对游戏开发感兴趣的开发者，还是对量化金融好奇的分析师，亦或是任何想直观理解数学如何解决实际问题的学习者，这都是一次绝佳的入门之旅。

2. 核心思路：微分方程的两种“解法”

在深入代码之前，我们必须统一思想：微分方程描述的是变化率。比如，速度是位置的变化率（导数），加速度是速度的变化率（二阶导数）。我们的核心任务，就是根据给定的变化率规则（微分方程本身）和初始状态（比如小球一开始的位置和速度），计算出物体在未来每个时刻的状态。

在计算机中，我们无法像解数学题那样得到一个完美的、连续的解析解。我们只能进行“数值求解”，即把连续的时间切成一片片很薄的时间片（时间步长dt），然后像做填空题一样，一片一片地推算出下一个时刻的状态。这里就有两种主流的“填空”思路，对应着我们两个实践场景。

2.1 显式欧拉法：游戏物理中的“快糙猛”

在游戏物理引擎中，如Unity内置的PhysX或自写的简单物理模拟，最常用的是显式欧拉法。它的思想直白得惊人：用当前时刻的状态，直接乘以变化率，再乘以时间步长，就得到了下一个时刻的状态。

公式看起来很简单：新位置 = 旧位置 + 速度 * dt，新速度 = 旧速度 + 加速度 * dt。

在Unity中，你几乎每天都在用它，只是你没察觉。在Update()或FixedUpdate()函数里，你写transform.position += velocity * Time.deltaTime;这就是最标准的显式欧拉积分，用于更新位置。

为什么游戏爱用它？

1. 计算简单，速度极快：一次乘法和一次加法，对需要每帧计算成千上万个物体的游戏来说，效率就是生命。
2. 易于实现和理解：代码直观，符合直觉，非常适合实时交互场景。

但它有个致命的“坑”：不稳定。特别是当系统刚度很大（比如一个很硬的弹簧）或者时间步长dt设置得不够小时，能量会越算越多，导致模拟爆炸（比如弹簧疯狂抖动直至飞出去）。游戏引擎通过一些技巧（如约束求解器、碰撞检测后的冲量修正）来弥补这个缺陷，但对于精度要求高的科学计算，这就力不从心了。

2.2 龙格-库塔法：金融模型中的“精细稳”

在金融、航天等对精度要求更高的领域，龙格-库塔法（尤其是四阶，简称RK4）是数值求解的黄金标准。它的思路比欧拉法“狡猾”得多：我不相信只用起点的一个斜率来预测整个步长内的变化。RK4会在时间步长内选取多个点（通常是起点、中点、终点）估算斜率，然后对这些斜率进行加权平均，得到一个更靠谱的“平均斜率”，再用这个平均斜率去更新状态。

为什么金融模型选它？

1. 精度高：对于光滑的系统，RK4的精度远高于欧拉法。在模拟资产价格这种对初始条件敏感（混沌性）的过程中，高精度能减少误差累积，让模型更可靠。
2. 稳定性好：能够容忍比显式欧拉法更大的时间步长，在长期模拟中不易发散。

代价是什么？计算量是显式欧拉法的四倍左右。但对于金融模型，我们通常不需要像游戏那样每秒模拟60帧，可能一天甚至一个月才走一步，计算量完全可接受。

注意：这里的选择不是绝对的。现代游戏物理引擎的底层约束求解器非常复杂，并非简单的显式欧拉。而一些简单的金融模型也可能用欧拉法。但对于入门理解，掌握这两种典型的、代表不同哲学的方法，至关重要。

3. Unity实战：构建一个微分方程可视化实验室

理论说再多，不如动手做。我们将在Unity里创建一个最小化的物理模拟场景，抛开复杂的引擎内部机制，亲手用代码实现微分方程的求解，并可视化每一步的过程。

3.1 环境与场景搭建

首先，打开Unity Hub，创建一个新的3D核心项目。在场景中，我们创建几个简单的物体：

1. 一个Cube：重命名为Floor，缩放其Scale为(10, 0.1, 10)，作为地面。
2. 一个Sphere：重命名为Ball，这就是我们模拟的主角。将其Y轴位置抬高，比如设为(0, 5, 0)。
3. 一个空物体：重命名为SimulationController，我们将把控制脚本挂在这里。

为Ball添加一个刚体组件Rigidbody，但关键一步来了：在Inspector面板中，找到Rigidbody组件，将Use Gravity（使用重力）取消勾选。因为我们要用自己的微分方程来模拟重力，而不是用Unity内置的物理引擎。

3.2 编写显式欧拉积分器脚本

在SimulationController上创建一个新的C#脚本，命名为ExplicitEulerSimulator。这个脚本将驱动小球的运动。

using UnityEngine; public class ExplicitEulerSimulator : MonoBehaviour { // 模拟的目标物体 public GameObject targetBall; // 重力加速度 (m/s^2)，负号表示方向向下 public float gravity = -9.81f; // 模拟的时间步长 (秒)。FixedUpdate默认是0.02s (50Hz)，这里我们显式控制。 public float simulationDeltaTime = 0.02f; // 当前速度 private Vector3 velocity = Vector3.zero; // 是否开始模拟 private bool isSimulating = false; void Start() { if (targetBall == null) { Debug.LogError("请指定要模拟的小球！"); return; } // 初始速度可以设为零，或者给一个初速度 velocity = Vector3.zero; isSimulating = true; } void FixedUpdate() { if (!isSimulating) return; // 1. 根据当前状态计算加速度（变化率）。 // 在这个自由落体例子中，加速度就是重力加速度，方向向下。 Vector3 acceleration = new Vector3(0, gravity, 0); // 2. 执行显式欧拉积分： // 新速度 = 旧速度 + 加速度 * dt velocity += acceleration * simulationDeltaTime; // 新位置 = 旧位置 + 速度 * dt targetBall.transform.position += velocity * simulationDeltaTime; // 3. 简单的地面碰撞检测（Y=0处） if (targetBall.transform.position.y <= 0.5f) // 假设小球半径0.5 { // 碰到地面，反转Y方向速度，并乘以一个弹性系数（小于1表示有能量损失） velocity.y = -velocity.y * 0.7f; // 防止小球钻入地下 Vector3 pos = targetBall.transform.position; pos.y = 0.5f; targetBall.transform.position = pos; } } // 提供一个方法，可以随时重置模拟 public void ResetSimulation(Vector3 startPosition) { if (targetBall != null) { targetBall.transform.position = startPosition; velocity = Vector3.zero; } } }

将脚本挂载到SimulationController上，并将场景中的Ball拖拽到脚本的Target Ball字段。运行游戏，你会看到小球做自由落体并弹跳，但运动轨迹可能看起来有点“硬”或者不自然，这就是显式欧拉和简单碰撞处理的特性。

3.3 引入阻尼振动与RK4对比

为了让例子更丰富，我们再模拟一个经典系统：弹簧振子。创建一个新的脚本SpringOscillator，这次我们同时实现显式欧拉和RK4，并对比效果。

using UnityEngine; public class SpringOscillator : MonoBehaviour { public enum IntegrationMethod { ExplicitEuler, RK4 } public IntegrationMethod method = IntegrationMethod.ExplicitEuler; // 弹簧参数 public float mass = 1.0f; // 质量 public float springConstant = 10.0f; // 弹簧系数 k，越大弹力越强 public float damping = 0.1f; // 阻尼系数 c，用于消耗能量 public Vector3 equilibriumPosition = Vector3.zero; // 平衡位置 private Vector3 position; // 当前位置 private Vector3 velocity; // 当前速度 private float simulationTime = 0f; public float dt = 0.02f; // 时间步长 void Start() { // 初始位置偏离平衡点，给予一个初始扰动 position = equilibriumPosition + new Vector3(2, 0, 0); velocity = Vector3.zero; simulationTime = 0f; } void FixedUpdate() { if (method == IntegrationMethod.ExplicitEuler) { UpdateExplicitEuler(); } else if (method == IntegrationMethod.RK4) { UpdateRK4(); } // 更新物体实际位置 transform.position = position; simulationTime += dt; } // 计算弹簧力：F = -k * (x - x0) - c * v （胡克定律 + 阻尼力） Vector3 CalculateAcceleration(Vector3 pos, Vector3 vel) { Vector3 displacement = pos - equilibriumPosition; Vector3 springForce = -springConstant * displacement; Vector3 dampingForce = -damping * vel; Vector3 totalForce = springForce + dampingForce; // a = F / m return totalForce / mass; } void UpdateExplicitEuler() { Vector3 acceleration = CalculateAcceleration(position, velocity); velocity += acceleration * dt; position += velocity * dt; } void UpdateRK4() { // RK4需要计算四个斜率(k1, k2, k3, k4) Vector3 pos = position; Vector3 vel = velocity; // k1: 起点斜率 Vector3 a1 = CalculateAcceleration(pos, vel); Vector3 k1_v = a1 * dt; // 速度的增量k1 Vector3 k1_x = vel * dt; // 位置的增量k1 // k2: 用k1预测的中点斜率 Vector3 a2 = CalculateAcceleration(pos + k1_x * 0.5f, vel + k1_v * 0.5f); Vector3 k2_v = a2 * dt; Vector3 k2_x = (vel + k1_v * 0.5f) * dt; // k3: 用k2预测的另一个中点斜率 Vector3 a3 = CalculateAcceleration(pos + k2_x * 0.5f, vel + k2_v * 0.5f); Vector3 k3_v = a3 * dt; Vector3 k3_x = (vel + k2_v * 0.5f) * dt; // k4: 用k3预测的终点斜率 Vector3 a4 = CalculateAcceleration(pos + k3_x, vel + k3_v); Vector3 k4_v = a4 * dt; Vector3 k4_x = (vel + k3_v) * dt; // 加权平均，得到最终增量 Vector3 deltaVelocity = (k1_v + 2f*k2_v + 2f*k3_v + k4_v) / 6f; Vector3 deltaPosition = (k1_x + 2f*k2_x + 2f*k3_x + k4_x) / 6f; // 更新状态 velocity += deltaVelocity; position += deltaPosition; } void OnGUI() { // 简单GUI切换方法 if (GUI.Button(new Rect(10, 10, 150, 30), "切换为显式欧拉")) method = IntegrationMethod.ExplicitEuler; if (GUI.Button(new Rect(10, 50, 150, 30), "切换为RK4")) method = IntegrationMethod.RK4; GUI.Label(new Rect(10, 90, 250, 30), $"当前方法: {method}"); GUI.Label(new Rect(10, 120, 250, 30), $"能量（近似）: {0.5f*mass*velocity.sqrMagnitude + 0.5f*springConstant*(position-equilibriumPosition).sqrMagnitude:F2}"); } }

将这个脚本挂载到一个新的小球上。运行后，你可以通过屏幕按钮切换两种积分方法。仔细观察：

• 显式欧拉：当弹簧系数k较大或阻尼c较小时，振动可能会逐渐放大（能量增加），导致模拟不稳定甚至崩溃。尝试把springConstant调到50，damping调到0.01看看。
• RK4：在相同参数下，振动通常更稳定，能量衰减更符合物理规律（阻尼消耗能量）。即使增大时间步长dt，RK4也比欧拉法更不容易发散。

这个对比实验直观地展示了不同数值方法在精度和稳定性上的巨大差异，这是选择算法的核心依据。

4. Python实战：用微分方程模拟金融资产价格

在Unity里我们模拟了确定性的物理系统（给定公式，结果唯一）。金融世界则充满了不确定性，资产价格的变化常被建模为随机微分方程。最经典的模型是几何布朗运动，它假设价格的对数收益率服从正态分布。其离散形式（欧拉-丸山格式）对我们来说已经足够入门。

4.1 环境准备与模型原理

我们使用Python，主要依赖numpy进行数值计算，matplotlib进行绘图。如果你还没有环境，建议使用Anaconda或直接pip install numpy matplotlib。

几何布朗运动的随机微分方程是：dS_t = μ * S_t * dt + σ * S_t * dW_t其中：

• S_t：时刻t的资产价格。
• μ：漂移率，代表资产的预期收益率。
• σ：波动率，代表资产价格的不确定性（风险）。
• dW_t：维纳过程（布朗运动）的增量，服从均值为0、方差为dt的正态分布。

这个方程的解（解析解）是已知的，但为了与我们之前的数值求解思想衔接，我们先用欧拉-丸山离散化来数值求解它，其离散形式为：S_{t+1} = S_t * (1 + μ * dt + σ * sqrt(dt) * Z)其中Z是标准正态分布随机数。

4.2 编写资产价格模拟器

创建一个新的Python文件，如financial_model.py。

import numpy as np import matplotlib.pyplot as plt import seaborn as sns sns.set_style("whitegrid") # 设置好看的绘图风格 def simulate_gbm_euler(S0, mu, sigma, T, dt, num_simulations=1): """ 使用欧拉-丸山方法模拟几何布朗运动（资产价格路径）。 参数: S0: 初始价格 mu: 年化漂移率（预期收益率） sigma: 年化波动率 T: 模拟总时长（年） dt: 时间步长（年） num_simulations: 模拟路径条数 返回: time_array: 时间序列 price_paths: 模拟的价格路径，形状为 (num_steps+1, num_simulations) """ num_steps = int(T / dt) # 总步数 time_array = np.linspace(0, T, num_steps + 1) # 时间轴 # 初始化价格路径矩阵，第一行都是初始价格S0 price_paths = np.zeros((num_steps + 1, num_simulations)) price_paths[0, :] = S0 # 生成随机数：标准正态分布 # 注意：这里使用向量化操作，一次生成所有需要的随机数，效率远高于循环 random_shocks = np.random.standard_normal((num_steps, num_simulations)) # 欧拉-丸山迭代 for t in range(1, num_steps + 1): # 关键公式：S_t = S_{t-1} * (1 + mu*dt + sigma*sqrt(dt)*Z) price_paths[t, :] = price_paths[t-1, :] * (1 + mu * dt + sigma * np.sqrt(dt) * random_shocks[t-1, :]) return time_array, price_paths # 模拟参数设置（以股票为例） S0 = 100.0 # 初始价格100元 mu = 0.08 # 预期年化收益率8% sigma = 0.20 # 年化波动率20%（典型股票水平） T = 2.0 # 模拟2年 dt = 1.0/252 # 时间步长：假设一年252个交易日，即模拟每个交易日 # 模拟5条可能的未来价格路径 time, paths = simulate_gbm_euler(S0, mu, sigma, T, dt, num_simulations=5) # 绘图 plt.figure(figsize=(12, 6)) for i in range(paths.shape[1]): plt.plot(time, paths[:, i], lw=1.5, alpha=0.7, label=f'Path {i+1}' if i < 3 else "") plt.title('Geometric Brownian Motion - Simulated Stock Price Paths (Euler-Maruyama)', fontsize=14) plt.xlabel('Time (Years)') plt.ylabel('Stock Price') plt.axhline(y=S0, color='gray', linestyle='--', alpha=0.5, label='Initial Price') plt.legend() plt.show()

运行这段代码，你会得到5条从同一起点出发，但走势各异的资产价格曲线。这就是随机微分方程的魅力：它描述了概率分布，而非单一确定路径。mu决定了曲线的整体向上趋势，而sigma决定了曲线的“毛刺”程度，即风险大小。

4.3 从模拟到应用：期权定价初探

模拟出价格路径有什么用？一个最经典的应用是期权定价。欧式看涨期权允许持有者在到期日以约定的行权价K购买资产。其到期收益为max(S_T - K, 0)。根据风险中性定价理论，期权的公平价格是其未来收益的期望值按无风险利率r折现到今天。

我们可以用刚才的模拟来近似计算这个期望值，这种方法称为蒙特卡洛模拟。

def monte_carlo_option_price(S0, K, T, r, sigma, num_simulations=100000): """ 使用蒙特卡洛模拟计算欧式看涨期权价格。 假设风险中性世界，即漂移率 mu = 无风险利率 r。 """ # 在风险中性测度下，漂移率等于无风险利率 mu_neutral = r dt = T / 100.0 # 将到期时间分为100步，精度足够 # 模拟大量路径在到期日T的价格 _, paths = simulate_gbm_euler(S0, mu_neutral, sigma, T, dt, num_simulations) # 取最后一天的价格 ST = paths[-1, :] # 计算每条路径的到期收益 payoffs = np.maximum(ST - K, 0) # 计算收益的期望值，并用无风险利率折现 option_price = np.exp(-r * T) * np.mean(payoffs) # 计算标准误，评估模拟精度 standard_error = np.std(payoffs) / np.sqrt(num_simulations) return option_price, standard_error # 期权参数 K = 105.0 # 行权价 r = 0.05 # 无风险利率5% num_sims = 50000 # 模拟次数 price_mc, se = monte_carlo_option_price(S0, K, T, r, sigma, num_sims) print(f"蒙特卡洛模拟的期权价格: {price_mc:.4f}") print(f"标准误: {se:.6f}") print(f"95% 置信区间: [{price_mc - 1.96*se:.4f}, {price_mc + 1.96*se:.4f}]") # （可选）与著名的Black-Scholes解析解公式对比 from scipy.stats import norm def black_scholes_call(S, K, T, r, sigma): d1 = (np.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma * np.sqrt(T)) d2 = d1 - sigma * np.sqrt(T) call_price = S * norm.cdf(d1) - K * np.exp(-r*T) * norm.cdf(d2) return call_price price_bs = black_scholes_call(S0, K, T, r, sigma) print(f"Black-Scholes 解析解价格: {price_bs:.4f}") print(f"差异: {abs(price_mc - price_bs):.6f}")

你会看到，通过模拟几万甚至几十万条价格路径，我们计算出的期权价格与Black-Scholes公式的解析解非常接近。这验证了我们模型的正确性，也展示了蒙特卡洛方法的强大之处：对于没有解析解的复杂衍生品，它几乎是唯一的通用定价工具。

5. 核心环节：参数、稳定性与性能的权衡

无论是游戏物理还是金融模型，当你真正动手实现时，都会遇到三个灵魂拷问：参数怎么设？模拟为什么炸了？跑得太慢了怎么办？

5.1 参数选择的艺术与科学

在Unity物理模拟中：

• 时间步长dt：这是最重要的参数。在FixedUpdate中，Unity默认Time.fixedDeltaTime为0.02秒（50Hz）。对于快速移动的物体或刚性碰撞，你可能需要更小的值（如0.01甚至0.005）来保证稳定性。但更小的dt意味着更高的CPU开销。经验法则：从默认值开始，如果出现物体穿透、剧烈抖动，就尝试减小dt。
• 迭代次数：对于复杂的约束（如多个关节连接），物理引擎内部需要迭代求解。在Project Settings -> Physics中增加Solver Iteration Count可以提高稳定性，但同样消耗性能。
• 质量、阻尼、弹簧常数：这些是模型参数。一个常见的坑是单位不统一。你设定的重力是-9.81，是米/秒^2。那么你物体的尺寸（Scale）最好接近真实世界（1单位=1米），质量（Mass）也应在合理范围（如一个人形角色70-100），否则会导致奇怪的物理行为。

在金融随机模拟中：

• 时间步长dt：对于几何布朗运动，离散化误差与dt成正比。通常，dt取为1/252（日度）或1/12（月度）对于大多数应用足够了。若要模拟高频交易，则需要更小的步长。
• 漂移率mu和波动率sigma：这些需要从历史数据中估计。mu的估计非常不准，且对长期预测影响巨大。在实践中，尤其是期权定价中，我们常使用风险中性测度，此时mu被替换为无风险利率r，从而规避了对未来预期收益的猜测。
• 模拟路径数num_simulations：蒙特卡洛模拟的精度与1/sqrt(N)成正比。要將标准误减半，你需要将模拟次数增加到4倍。通常，10万次模拟对于普通期权定价已经能提供不错的结果（标准误很小）。需要在精度和计算时间之间权衡。

5.2 数值稳定性：为什么我的模拟爆炸了？

这是数值计算永恒的课题。

1. 显式方法的条件稳定性：如前所述，显式欧拉法是有条件稳定的。对于方程dx/dt = -k*x（类似衰减），要保证稳定，需要k*dt < 2。对于弹簧振子，稳定条件与sqrt(k/m)*dt有关。对策：减小时间步长dt是首选。如果无法减小，则需换用隐式方法（如后向欧拉），这类方法无条件稳定但计算更复杂，Unity的某些物理求解器内部会使用。

2. 随机模拟中的负价格：在几何布朗运动的欧拉离散化中，如果波动率sigma很大，而时间步长dt不够小，随机冲击sigma * sqrt(dt) * Z可能小于-1，导致1 + mu*dt + ...为负，从而产生负的价格，这在金融上是不可能的。对策：使用对数欧拉方法，即对价格取对数进行离散化，或者直接使用解析解的离散形式S_{t+1} = S_t * exp( (mu - 0.5*sigma^2)*dt + sigma*sqrt(dt)*Z )，这能保证价格永远为正。

def simulate_gbm_analytic(S0, mu, sigma, T, dt, num_simulations): """使用几何布朗运动解析解的离散形式模拟，保证价格为正。""" num_steps = int(T / dt) price_paths = np.zeros((num_steps + 1, num_simulations)) price_paths[0, :] = S0 # 一次生成所有随机数 Z = np.random.standard_normal((num_steps, num_simulations)) for t in range(1, num_steps + 1): # 解析解离散形式 price_paths[t, :] = price_paths[t-1, :] * np.exp( (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z[t-1, :] ) return price_paths

5.3 性能优化技巧

Unity端：

• 减少FixedUpdate负载：只在必要时进行复杂的自定义物理计算。对于大量静态或简单运动的物体，使用Unity原生物理并做好碰撞层管理。
• 对象池：对于频繁生成和销毁的抛射物等，使用对象池复用GameObject，避免昂贵的实例化和垃圾回收。
• 批处理与并行：对于大量独立物体的模拟（比如一群鸟），如果自定义逻辑，可以考虑使用Job System和Burst Compiler进行C# Job多线程并行计算，性能提升显著。

Python/金融端：

• 向量化操作：永远避免在Python中使用显式循环遍历数组。像之前代码中用np.random.standard_normal一次性生成所有随机数，用数组运算更新所有路径，就是向量化。对比循环版本，速度可能有百倍差距。
• 使用numba加速：对于无法向量化的复杂逻辑，可以使用@numba.jit装饰器将函数编译为机器码。对于蒙特卡洛模拟，这通常能带来数十倍的加速。
• 并行计算：如果模拟路径之间完全独立，可以轻松并行。使用concurrent.futures或joblib库，或者numba的并行模式。

import numba import numpy as np @numba.jit(nopython=True, parallel=True) # 启用并行 def simulate_gbm_numba(S0, mu, sigma, T, dt, num_simulations): num_steps = int(T / dt) price_paths = np.zeros((num_steps + 1, num_simulations)) price_paths[0, :] = S0 # Numba 的并行循环 for i in numba.prange(num_simulations): for t in range(1, num_steps + 1): Z = np.random.normal() # Numba内部优化 price_paths[t, i] = price_paths[t-1, i] * np.exp( (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z ) return price_paths # 首次运行会有编译开销，后续运行极快

6. 常见问题与调试心得

在实际操作中，你肯定会遇到各种奇怪的问题。这里记录一些我踩过的坑和解决方法。

Unity物理模拟部分：

1. 小球穿墙而过

   • 问题：速度太快，一帧内移动距离超过了碰撞体厚度。
   • 解决：Unity的Continuous碰撞检测模式就是为此设计的。在Rigidbody组件中，将Collision Detection从Discrete改为Continuous或Continuous Dynamic。或者，在自己实现的简单模拟中，进行更精细的“连续碰撞检测”，比如计算本帧移动的射线，而不是只检查终点。

2. 模拟抖动（Jitter）

   • 问题：物体，尤其是堆叠的物体，不停微幅震动。
   • 解决：增加物理引擎的求解器迭代次数（Solver Iteration Count）。适当增加物体的质量（Mass），减少极端质量比（如一个质量1的物体去撞质量10000的物体）。检查是否有多余的力在持续作用。

3. 自定义物理与Unity物理混合使用导致鬼畜

   • 问题：你用自己的代码控制位置，但物体上又有Rigidbody，Unity物理引擎也在计算。
   • 解决：二选一。要么彻底禁用Unity物理（Rigidbody设置为Is Kinematic，自己处理所有碰撞），要么完全交给Unity物理，通过AddForce等方式施加力。混合模式需要极其小心地协调，通常不建议。

Python金融模拟部分：

1. 模拟结果每次运行都不一样

   • 问题：这是正常的！因为使用了随机数。为了结果可复现，需要在模拟前设置随机数种子。
   • 解决：在代码开头添加np.random.seed(42)# 42是任意选的种子。这样每次运行都会生成相同的随机序列。

2. 蒙特卡洛模拟价格与解析解差距较大

   • 检查1：模拟次数是否足够？计算标准误，看置信区间是否包含解析解。
   • 检查2：离散化方法是否正确？对于几何布朗运动，优先使用基于解析解的离散形式（exp那个），避免欧拉离散可能带来的偏差。
   • 检查3：参数单位是否一致？mu,sigma,r,T必须使用相同的时间单位（通常都是年化）。如果T是2年，那么mu=0.08代表年化8%。

3. 代码运行太慢

   • 定位：用%%timeit（Jupyter）或time模块测量函数各部分耗时。瓶颈通常在循环和随机数生成。
   • 优化：如前所述，矢量化是第一步。对于超大规模模拟（如百万路径），考虑使用numba或Cython，甚至调用C++库。

一个通用的调试心法：从简单到复杂，逐步验证。

1. 先去掉所有随机性（设sigma=0），看确定性部分（mu的影响）是否正确。
2. 再测试随机性，可以固定随机数种子，对比多次运行结果。
3. 对于物理模拟，先在一个静止状态下开始，然后施加一个恒力，看运动是否符合F=ma的预期。
4. 大量使用可视化。在Unity里画Debug射线、向量；在Python里把中间变量（如每次的收益率、随机冲击）都画出来看看分布是否符合预期。图形比数字更能直观地暴露问题。

从Unity中一个球的下落，到Python里一条资产价格曲线的生成，微分方程就像一条隐藏的丝线，串起了虚拟与现实中的动态系统。这次实践的核心收获，不是记住了欧拉或龙格-库塔的公式，而是建立起一种“数值思维”：面对一个描述变化率的方程，我知道如何将它拆解成计算机能执行的微小步骤，并清醒地意识到每种拆解方法带来的精度、稳定性和性能上的取舍。这种思维，是你在游戏里设计一个新颖的物理交互（比如一个基于流体的魔法效果），或者在金融领域构建一个更复杂的随机波动率模型时，真正赖以创新的基础。工具和库会变，但这种将连续世界离散化的能力，会一直是你工具箱里最趁手的一件。