26、伊藤随机微积分：随机微分方程及其数值解法

伊藤随机微积分：随机微分方程及其数值解法

1. 随机微分方程概述

随机微分方程在描述具有随机因素的动态系统中具有重要作用。我们将重点关注朗之万方程，它是描述布朗粒子速度的模型，类似于在常微分方程研究中简谐运动的作用，能以不过于复杂的方式阐释随机微分方程的诸多方面。

1.1 利用乘积法则求解的示例

1.1.1 示例一

考虑随机微分方程：
[dX(t) = [t + B^2(t)] dt + 2tB(t) dB(t), \quad X(0) = X_0]
通过巧妙运用乘积法则，我们发现：
[dX(t) = B^2(t) dt + t[2B(t) dB(t) + dt] = B^2(t) dt + t d[B^2(t)] = d[tB^2(t)]]
对两边进行积分，可得该方程的解为：
[X(t) = tB^2(t) + X_0]

1.1.2 示例二

对于随机微分方程：
[dX(t) = \frac{b - X(t)}{1 - t} dt + dB(t), \quad 0 \leq t < 1, \quad X(0) = X_0]
首先将其改写为：
[\frac{d[b - X(t)]}{1 - t} + \frac{b - X(t)}{(1 - t)^2} dt = -\frac{dB(t)}{1 - t}]
逆向运用乘积法则可得：
[d\left[\frac{b - X(t)}{1 - t}\right] = -\frac{dB(t)}{1 - t}]
从 0 到 t 对两边积分：
[\frac{b - X(t)}{1 - t} = b - X(0) - \int_0^t \frac{dB(\eta)}{1 - \eta}]
求解(X(t))，最终结果为：
[X(t) = b - b - X(0) + (1 - t) \int_0^t \frac{dB(\eta)}{1 - \eta}]
在这种情况下，若需要数值解，由于积分难以化简，需应用数值求积法。

1.2 朗之万方程的求解

朗之万方程为：
[dX(t) = cX(t) dt + \sigma dB(t), \quad X(0) = X_0]
其解为：
[X(t) = X_0 + c \int_0^t X(\eta) d\eta + \sigma \int_0^t dB(\eta)]
为求解该方程，我们引入函数(f(t, x) = e^{-ct}x)，利用伊藤引理可得：
[f[t, X(t)] - X(0) = \int_0^t \left{f_t[\eta, X(\eta)] + cX(\eta)f_x[\eta, X(\eta)] + \frac{1}{2}\sigma^2f_{xx}[\eta, X(\eta)]\right} d\eta + \int_0^t \sigma f_x[\eta, X(\eta)] dB(\eta)]
将(f(t, x))代入上式，经过化简可得：
[e^{-ct}X(t) - X_0 = \sigma \int_0^t e^{-c\eta} dB(\eta)]
求解(X(t))，得到：
[X(t) = X(0)e^{ct} + \sigma e^{ct} \int_0^t e^{-c\eta} dB(\eta)]
当(X_0)为常数时，(X(t))被称为奥恩斯坦 - 乌伦贝克过程。

另一种推导方法是将朗之万方程乘以积分因子(e^{-ct})，得到：
[e^{-ct} dX(t) - ce^{-ct}X(t) dt = \sigma e^{-ct} dB(t)]
逆向运用乘积法则可得：
[d\left[e^{-ct}X(t)\right] = \sigma e^{-ct} dB(t)]
对两边积分，同样可得到上述结果。

1.3 精确随机微分方程示例

考虑随机微分方程：
[X(t) = X(0) + c \int_0^t X(s) ds + \sigma \int_0^t X(s) dB(s), \quad c, \sigma > 0]
若(X(t) = f[t, B(t)])，根据伊藤引理可得：
[X(t) = X(0) + \int_0^t \left{f_t[s, B(s)] + \frac{1}{2}f_{xx}[s, B(s)]\right} ds + \int_0^t f_x[s, B(s)] dB(s)]
对比两个方程，我们得到：
[\begin{cases}
cf(t, x) = f_t(t, x) + \frac{1}{2}f_{xx}(t, x) \
\sigma f(t, x) = f_x(t, x)
\end{cases}]
由第二个方程可得(f_{xx}(t, x) = \sigma f_x(t, x) = \sigma^2 f(t, x))，代入第一个方程可得：
[(c - \frac{1}{2}\sigma^2) f(t, x) = f_t(t, x)]
使用分离变量法求解上述方程组，可得：
[f(t, x) = f(0, 0) \exp\left[(c - \frac{1}{2}\sigma^2) t + \sigma x\right]]
即：
[X(t) = f[t, B(t)] = X(0) \exp\left[(c - \frac{1}{2}\sigma^2) t + \sigma B(t)\right]]
这种解被称为几何布朗运动。

1.4 齐次线性方程示例

对于齐次线性随机微分方程：
[dX(t) = c_1(t)X(t) dt + \sigma_1(t)X(t) dB(t)]
引入函数(f(t, x) = \ln(x))，根据伊藤引理可得：
[df = d[\ln(X)] = \left[c_1(t) - \frac{1}{2}\sigma_1^2(t)\right] dt + \sigma_1(t) dB(t)]
对两边积分并取指数，得到该方程的解为：
[X(t) = X(0) \exp\left[\int_0^t \left[c_1(\eta) - \frac{1}{2}\sigma_1^2(\eta)\right] d\eta + \int_0^t \sigma_1(\eta) dB(\eta)\right]]

1.5 一般情况示例

考虑齐次线性随机微分方程：
[dX(t) = [c_1(t)X(t) + c_2(t)] dt + [\sigma_1(t)X(t) + \sigma_2(t)] dB(t)]
我们先考虑齐次线性随机微分方程：
[dY(t) = c_1(t)Y(t) dt + \sigma_1(t)Y(t) dB(t), \quad Y(0) = 1]
由前面的示例可知，其解为：
[Y(t) = \exp\left{\int_0^t \left[c_1(\eta) - \frac{1}{2}\sigma_1^2(\eta)\right] d\eta + \int_0^t \sigma_1(\eta) dB(\eta)\right}]
引入两个随机变量(X_1 = \frac{1}{Y})和(X_2 = X)，利用伊藤引理可得：
[dX_1 = \left[\sigma_1^2(t) - c_1(t)\right] X_1(t) dt - \sigma_1(t)X_1(t) dB(t)]
再根据相关公式可得：
[d(X_1X_2) = [c_2(t) - \sigma_1(t)\sigma_2(t)] X_1(t) dt + \sigma_2(t)X_1(t) dB(t)]
对两边积分，最终得到：
[X(t) = Y(t) \left{X(0) + \int_0^t \frac{[c_2(\eta) - \sigma_1(\eta)\sigma_2(\eta)]}{Y(\eta)} d\eta + \int_0^t \frac{\sigma_2(\eta) dB(\eta)}{Y(\eta)}\right}]

1.6 随机 Verhulst 方程示例

随机 Verhulst 方程为：
[dX(t) = aX(t)[M - X(t)] dt + bX(t) dB(t), \quad X(0) = X_0]
引入(\Phi(t) = \frac{1}{X(t)})，根据伊藤引理可得：
[d\Phi(t) = -\Phi(t)[(aM - b^2) dt + b dB(t)] + a dt, \quad \Phi(0) = \frac{1}{X_0}]
利用前面示例的结果求解该方程，最终可得：
[X(t) = \frac{X_0 \exp[\xi(t)]}{1 + aX_0 \int_0^t \exp[\xi(\eta)] d\eta}]
其中(\xi(t) = (aM - \frac{b^2}{2})t + bB(t))。

1.7 练习题

以下是一些练习题，用于巩固所学知识：
1. 求解随机微分方程(dX(t) = \frac{1}{2}e^{t/2}B(t) dt + e^{t/2} dB(t), \quad X(0) = X_0)，通过逆向运用乘积法则。
2. 求解随机微分方程(dX(t) = e^{2t}[1 + 2B^2(t)] dt + 2e^{2t}B(t) dB(t), \quad X(0) = X_0)，提示：重写微分方程为(dX(t) = e^{2t}[2B(t) dB(t) + dt] + (2e^{2t} dt)B^2(t))。
3. 求解随机微分方程(dX(t) = [1 + B(t)] dt + [t + 2B(t)] dB(t), \quad X(0) = X_0)，提示：重写微分方程为(dX(t) = 2B(t) dB(t) + dt + B(t) dt + t dB(t))。
4. 求解随机微分方程(dX(t) = [3t^2 + B(t)] dt + t dB(t), \quad X(0) = X_0)，提示：重写微分方程为(dX(t) = 3t^2 dt + [B(t) dt + t dB(t)])。
5. 求解随机微分方程(dX(t) = B^2(t) dt + 2tB(t) dB(t), \quad X(0) = X_0)，提示：重写微分方程为(dX(t) = t[2B(t) dB(t) + dt] + B^2(t) dt - t dt)。
6. 求随机微分方程(dX(t) = [\beta - \alpha X(t)] dt + \sigma dB(t), \quad X(0) = X_0)的积分因子和解，其中(B(t))是布朗运动，(\alpha)，(\beta)和(\sigma)是常数。
7. 求随机微分方程(dX(t) = [1 + 2X(t)] dt + e^{2t} dB(t), \quad X(0) = X_0)的积分因子和解，其中(B(t))是布朗运动。
8. 求随机微分方程(dQ(t) + \frac{Q(t)}{RC} dt = \frac{V(t)}{R} dt + \frac{\alpha(t)}{R} dB(t), \quad Q(0) = Q_0)的积分因子和解，其中(R)和(C)是正实数常数，(B(t))是布朗运动。
9. 求随机微分方程(dX(t) = 2tX(t) dt + e^{-t} dt + dB(t), \quad t \in [0, 1], \quad X(0) = X_0)的积分因子和解，其中(B(t))是布朗运动。
10. 求随机微分方程(dX(t) = [4X(t) - 1] dt + 2 dB(t), \quad X(0) = X_0)的积分因子和解，其中(B(t))是布朗运动。
11. 求随机微分方程(dX(t) = [2 - X(t)] dt + e^{-t}B(t) dB(t), \quad X(0) = X_0)的积分因子和解，其中(B(t))是布朗运动。
12. 求随机微分方程(dX(t) = [1 + X(t)] dt + e^{t}B(t) dB(t), \quad X(0) = X_0)的积分因子和解，其中(B(t))是布朗运动。
13. 求随机微分方程(dX(t) = \left[\frac{1}{2}X(t) + 1\right] dt + e^{t} \cos[B(t)] dB(t), \quad X(0) = X_0)的积分因子和解，其中(B(t))是布朗运动。
14. 求随机微分方程(dX(t) = \left[t + \frac{1}{2}X(t)\right] dt + e^{t} \sin[B(t)] dB(t), \quad X(0) = X_0)的积分因子和解，其中(B(t))是布朗运动。
15. 仿照示例 8.5.2，求解精确随机微分方程(dX(t) = e^{t} [1 + B^2(t)] dt + [1 + 2e^{t}B(t)] dB(t), \quad X(0) = X_0)，步骤如下：
- 证明(f_t + \frac{1}{2}f_{xx} = e^{t}(1 + x^2))，且(f_x = 1 + 2e^{t}x)。
- 证明(f(t, x) = x + e^{t}x^2 + g(t))。
- 证明(g(t) = X_0)，且(X(t) = B(t) + e^{t}B^2(t) + X_0)。
16. 仿照示例 8.5.2，求解精确随机微分方程(dX(t) = \left{2tB^2(t) + 3t^2 [1 + B(t)]\right} dt + [1 + 3t^2B^2(t)] dB, \quad X(0) = X_0)，步骤如下：
- 证明(f_t + \frac{1}{2}f_{xx} = 2tx^3 + 3t^2(1 + x))，且(f_x = 3t^2x^2 + 1)。
- 证明(f(t, x) = t^2x^3 + x + g(t))。
- 证明(g’(t) = 3t^2)。
- 证明(X(t) = t^2[B^3(t) + t] + B(t) + X_0)。
17. 利用公式(X(t) = X(0) \exp\left[\int_0^t \left[c_1(\eta) - \frac{1}{2}\sigma_1^2(\eta)\right] d\eta + \int_0^t \sigma_1(\eta) dB(\eta)\right])，求解随机微分方程(dX(t) = t^2X(t) dt + tX(t) dB(t), \quad X(0) = X_0)。
18. 利用上述公式，求解随机微分方程(dX(t) = \cos(t)X(t) dt + \sin(t)X(t) dB(t), \quad X(0) = X_0)。
19. 利用上述公式，求解随机微分方程(dX(t) = \ln(t + 1)X(t) dt + \sqrt{\ln(t + 1)} X(t) dB(t), \quad X(0) = X_0)。
20. 利用上述公式，求解随机微分方程(dX(t) = \ln(t + 1)X(t) dt + tX(t) dB(t), \quad X(0) = X_0)。
21. 仿照示例 8.5.5，求解随机微分方程(dX(t) = [aX^n(t) + bX(t)] dt + cX(t) dB(t), \quad X(0) = X_0)，其中(n > 1)，步骤如下：
- 设(\Phi(t) = X^{1 - n}(t))，利用伊藤引理证明(d\Phi(t) = (1 - n)\Phi(t) \left[(b - \frac{1}{2}nc^2) dt + c dB(t)\right] + (1 - n)a dt)。
- 设(c_1(t) = (1 - n)b - \frac{n(1 - n)c^2}{2})，(c_2(t) = (1 - n)a)，(\sigma_1(t) = (1 - n)c)，(\sigma_2(t) = 0)，证明(\frac{\exp[(n - 1)\xi(t)]}{X^{n - 1}(t)} - \frac{1}{X_0^{n - 1}} = (1 - n)a \int_0^t \exp[(n - 1)\xi(\eta)] d\eta)，其中(\xi(t) = (b - \frac{c^2}{2})t + cB(t))。
22. 仿照示例 8.5.5，求解随机 Ginzburg - Landau 方程(dX(t) = \left[a e^{cX(t)} + b\right] dt + \sigma dB(t), \quad X(0) = X_0)，步骤如下：
- 设(\Phi(t) = \exp[-cX(t)])，利用伊藤引理证明(d\Phi(t) = -\left(bc - \frac{1}{2}\sigma^2c^2\right) \Phi(t) dt - \sigma c\Phi(t) dB(t) - ac dt)。
- 设(c_1(t) = \frac{\sigma^2c^2}{2} - bc)，(c_2(t) = -ac)，(\sigma_1(t) = -\sigma c)，(\sigma_2(t) = 0)，证明(X(t) = X_0 + bt + \sigma B(t) - \frac{1}{c} \ln\left{1 - ac \int_0^t \exp[cX_0 + bc\xi + \sigma cB(\xi)] d\xi\right})。
23. 仿照示例 8.5.5，求解随机微分方程(dX(t) = \left{[1 + X(t)][1 + X^2(t)]\right} dt + [1 + X^2(t)] dB(t), \quad X(0) = X_0)，步骤如下：
- 设(\Phi(t) = \tan^{-1}[X(t)])，利用伊藤引理证明(d\Phi(t) = dt + dB(t))。
- 求解步骤 1 中的随机微分方程，证明(X(t) = \tan[\tan^{-1}(X_0) + t + B(t)])。

1.8 求解流程总结

以下是利用乘积法则求解随机微分方程的一般流程：

graph TD;
    A[给定随机微分方程] --> B[尝试重写方程以运用乘积法则];
    B --> C{是否可运用乘积法则};
    C -- 是 --> D[逆向运用乘积法则化简方程];
    D --> E[对化简后的方程两边积分];
    E --> F[求解积分得到方程的解];
    C -- 否 --> G[考虑其他方法求解];

通过以上示例和练习题，我们可以看到利用乘积法则、伊藤引理等方法可以求解多种类型的随机微分方程。在实际应用中，需要根据具体方程的形式选择合适的方法进行求解。

2. 随机微分方程的数值解法

在实际应用中，很多随机微分方程无法通过解析方法求解，因此需要使用数值方法来近似求解。下面将介绍几种常见的数值解法。

2.1 欧拉 - 丸山（Euler - Marugama）方法

考虑随机微分方程：
[dX(t) = a[X(t), t] dt + b[X(t), t] dB(t)]
在区间(t_0 \leq t \leq T)上，初始值为(X(t_0) = X_0)。

我们引入网格(t_0 < t_1 < t_2 < \cdots < t_n < \cdots < t_N = T)，为了简化，假设所有时间增量相同，即(\Delta t = t_{n + 1} - t_n)，且(0 < \Delta t < 1)。

则有：
[X_{n + 1} = X_n + \int_{t_n}^{t_{n + 1}} a[X(\eta), \eta] d\eta + \int_{t_n}^{t_{n + 1}} b[X(\eta), \eta] dB(\eta)]
对上述积分进行最粗略的近似：
[\int_{t_n}^{t_{n + 1}} a[X(\eta), \eta] d\eta \approx a[X(t_n), t_n] \Delta t_n]
[\int_{t_n}^{t_{n + 1}} b[X(\eta), \eta] dB(\eta) \approx b[X(t_n), t_n] \Delta B_n]
其中(\Delta B_n = B(t_{n + 1}) - B(t_n))。

将这些近似代入(X_{n + 1})的表达式，得到欧拉 - 丸山近似：
[X_{n + 1} = X_n + a(t_n, X_n) (t_{n + 1} - t_n) + b(t_n, X_n) (B_{t_{n + 1}} - B_{t_n})]
对于(n = 0, 1, 2, \cdots, N - 1)，初始值为(X_0)。

当(b = 0)时，该随机迭代方案退化为常微分方程的常规欧拉方案。当(b \neq 0)时，由于(\Delta B_n)是独立的高斯随机变量，其均值(E(\Delta B_n) = 0)，方差(E[(\Delta B_n)^2] = \Delta t)。可以使用 MATLAB 函数 randn 来生成(\Delta B_n)。

2.2 收敛阶

在使用任何数值方案时，收敛速度是一个重要的考虑因素。在数值模拟过程中，在时间(t)处，精确解(X(t))和数值近似(Y(t))之间会存在差异(e(t) = X(t) - Y(t))，这也是一个随机变量。

如果对于任意时间(t)，当时间步长(\Delta t)足够小时，(E(|e(t)|) = O[(\Delta t)^m])，则称随机微分方程方案以阶(m)强收敛。可以证明，欧拉 - 丸山方法的强收敛阶为(\frac{1}{2})。

2.3 米尔斯坦（Milstein）方法

为了构造一个强收敛阶为 1 的近似方法，我们回到(X_{n + 1})的表达式：
[X_{n + 1} - X_n = \int_{t_n}^{t_{n + 1}} \left[a[X_n(\eta), \eta] + \int_{t_n}^{\eta} \left(a a_x + \frac{1}{2}b^2 a_{xx}\right) d\xi + \int_{t_n}^{\eta} b a_x dB(\xi)\right] d\eta + \int_{t_n}^{t_{n + 1}} \left[b[X_n(\eta), \eta] + \int_{t_n}^{\eta} \left(a b_x + \frac{1}{2}b^2 b_{xx}\right) d\xi + \int_{t_n}^{\eta} b b_x dB(\xi)\right] d\eta]
经过一系列化简和近似，忽略高阶项后，得到：
[R_n \approx b[X(t_n), t_n] b_x[X(t_n), t_n] \left[(\Delta B_n)^2 - \Delta t\right]]
最终得到米尔斯坦方法的表达式：
[X_{n + 1} = X_n + a(X_n, t_n) \Delta t_n + b(X_n, t_n) \Delta B_n + \frac{1}{2}b(X_n, t_n) \frac{\partial b(X_n, t_n)}{\partial x} \left[(\Delta B_n)^2 - \Delta t\right]]

2.4 示例：线性随机微分方程的数值解

考虑由线性随机微分方程定义的伊藤过程(X(t))：
[dX(t) = aX(t) dt + bX(t) dB(t)]
在区间(t \in [0, T])上，其精确解为：
[X(t) = X_0 \exp\left[\left(a - \frac{b^2}{2}\right) t + bB(t)\right]]

使用欧拉 - 丸山方法和米尔斯坦方法对该随机微分方程进行数值求解，并与精确解进行比较。不同时间步长下的数值解如图所示（此处可参考原文中的图 8.6.1）。

2.5 样本均值和标准差

虽然绘制各种实现可以大致了解随机过程对解的影响，但两个更有用的参数是在时间(t_n)处的样本均值和标准差：
- 样本均值：(\overline{X}(t_n) = \frac{1}{J} \sum_{j = 1}^{J} X_j(t_n))
- 样本标准差：(\sigma^2(t_n) = \frac{1}{J - 1} \sum_{j = 1}^{J} [X_j(t_n) - \overline{X}(t_n)]^2)
其中(J)是实现的数量，(X_j(t_n))是第(j)次实现中随机变量在时间(t_n)的值。

在许多物理问题中，“噪声”是随机过程的来源，我们假设服从正态分布(N(\mu, \sigma^2))，其中(\mu)和(\sigma)分别是总体均值和标准差。可以使用样本统计量来确定双侧置信区间：
[\left[\overline{X}(t_n) - \tau_{student} \frac{\sigma(t_n)}{\sqrt{J}}, \overline{X}(t_n) + \tau_{student} \frac{\sigma(t_n)}{\sqrt{J}}\right]]
基于自由度为(J - 1)的学生(-\tau)分布。

2.6 项目示例

2.6.1 RL 电气电路与噪声

考虑一个由电阻(R)和电感(L)组成的简单电气电路，由电压源(v(t))驱动，电流(I(t))满足一阶常微分方程：
[L \frac{dI}{dt} + RI = v(t), \quad I(0) = I_0]
其确定性解为：
[I(t) = I_0 e^{-\frac{Rt}{L}} + \frac{1}{L} \int_0^t \exp\left[\frac{R}{L} (\tau - t)\right] v(\tau) d\tau]

当考虑随机因素时，方程变为：
[dI = \frac{1}{L} [v(t) - RI(t)] dt - \frac{\alpha}{L} I(t) dB_1(t) + \frac{\beta}{L} dB_2(t), \quad I(0) = I_0]
其中(B_1(t))和(B_2(t))是两个独立的白噪声过程，(\alpha)，(\beta)是非负常数。

操作步骤如下：
1. 使用经典方法求解确定性方程，得到上述确定性解。
2. 使用 MATLAB 编写脚本，对随机微分方程进行数值积分，绘制多个实现下电流(I(t))随无量纲时间(\frac{Rt}{L})的变化曲线。
3. 计算并绘制解的均值和标准差随无量纲时间的变化曲线，并与确定性解进行比较。

2.6.2 具有布朗运动强迫的松弛振荡器

FitzHugh - Nagamo 模型描述了可激发系统，如神经元。修改后的模型方程为：
[\begin{cases}
dx = -x(x^2 - a^2) dt - y dt + \sigma dB_1(t) \
dy = \frac{(x - my) dt}{\tau} + \sigma dB_2(t)
\end{cases}]
其中(a)，(m)，(\sigma)和(\tau)是参数。

操作步骤如下：
1. 编写 MATLAB 脚本，对修改后的 FitzHugh - Nagamo 模型进行数值积分，针对不同的(\sigma)值进行多次模拟。
2. 计算(E[x(t)])随时间(t)的变化，分析布朗运动强迫的影响。

2.6.3 随机阻尼谐振子

随机阻尼谐振子由以下随机微分方程描述：
[\begin{cases}
dv(t) = -\gamma v(t) dt - k^2 x(t) dt - \alpha x(t) dB(t) \
dx(t) = v(t) dt
\end{cases}]
其中(k)，(\alpha)和(\gamma)是实常数。

可以使用欧拉方法或 Heun 方法进行数值求解。Heun 方法的步骤如下：
1. 首先进行欧拉式时间步长估计：
- (x^ = x_i + v_i \Delta t)
- (v^ = v_i - \gamma v_i \Delta t - k^2 x_i \Delta t - \alpha x_i \Delta B_i)
2. 然后使用估计值计算(x_{i + 1})和(v_{i + 1})：
- (x_{i + 1} = x_i + \frac{1}{2}(v_i + v^ ) \Delta t)
- (v_{i + 1} = v_i - \frac{1}{2} \gamma (v_i + v^ ) \Delta t - \frac{1}{2} k^2 (x_i + x^*) \Delta t - \alpha x_i \Delta B_i)

操作步骤如下：
1. 编写 MATLAB 脚本，使用欧拉方法和 Heun 方法对随机谐振子进行数值积分，设置参数(10\alpha = 4\gamma = k = 1)，初始条件(x(0) = v(0) = 0)。
2. 进行多次模拟，计算(E[x(t)])随时间(t)的变化，分析不同时间步长对解的精度的影响。

2.6.4 平均首次通过时间

随机微分方程：
[dX(t) = [X(t) - X^3(t)] dt + \frac{4}{X^2(t) + 1} dB(t)]
描述了在双势阱(V(x) = \frac{x^4}{4} - \frac{x^2}{2})中粒子的运动。

操作步骤如下：
1. 编写 MATLAB 代码，计算(X(t))随时间(t)的变化。
2. 进行(N)次实现，记录每次实现中粒子从初始位置(X(0) = -1)到达局部最大值(X(t) = 0)所需的时间。
3. 计算这些时间的平均值，并将结果绘制为时间步长平方根(\sqrt{h})的函数。

2.6.5 公司破产问题

随机微分方程：
[dX(t) = [\mu X(t) - iD] dt + \sigma X(t) dB(t), \quad 0 < t < T, \quad X(0) = X_0]
描述了公司财富(X(t))随时间(t)的演变。其中(\mu)和(\sigma)分别表示公司财富的确定性和随机演变，(X_0)是公司的初始财富，(iD)是向金融机构（银行）支付的款项。

操作步骤如下：
1. 编写 MATLAB 代码，模拟公司在给定(D)，(i)和(X_0)，以及不同(\sigma)值下的财富演变。
2. 在模拟过程中，记录公司破产（即(X(\tau) = 0)）的次数(n)，计算破产概率(P[X(\tau) = 0] = \frac{n}{N})。
3. 计算不同规模公司（小、中、大）的破产概率随利率的变化，分析平均破产时间(\tau)随利率的变化。

2.7 数值求解流程总结

以下是使用欧拉 - 丸山方法和米尔斯坦方法求解随机微分方程的一般流程：

graph TD;
    A[给定随机微分方程] --> B[确定网格和时间步长];
    B --> C[初始化初始值 X_0];
    C --> D[循环计算 X_{n + 1}];
    D --> E{选择方法};
    E -- 欧拉 - 丸山方法 --> F[使用欧拉 - 丸山公式计算 X_{n + 1}];
    E -- 米尔斯坦方法 --> G[使用米尔斯坦公式计算 X_{n + 1}];
    F --> H[更新 n = n + 1];
    G --> H;
    H --> I{是否达到最终时间 T};
    I -- 否 --> D;
    I -- 是 --> J[输出数值解];

通过以上介绍，我们了解了随机微分方程的数值解法，包括欧拉 - 丸山方法和米尔斯坦方法，以及它们在不同实际问题中的应用。在实际应用中，需要根据具体问题选择合适的数值方法，并注意时间步长对解的精度的影响。