2025-11-17T00:52:13.221997

On the random generation of Butcher trees

Huang, Privault
The main goal of this paper is to provide an algorithm for the random sampling of Butcher trees and the probabilistic numerical solution of ordinary differential equations (ODEs). This approach complements and simplifies a recent approach to the probabilistic representation of ODE solutions, by removing the need to generate random branching times. The random sampling of trees is compared to the finite order truncation of Butcher series in numerical experiments.
academic

On the random generation of Butcher trees

基本信息

  • 论文ID: 2404.05969
  • 标题: On the random generation of Butcher trees
  • 作者: Qiao Huang (Southeast University), Nicolas Privault (Nanyang Technological University)
  • 分类: math.CA (Classical Analysis), math.PR (Probability)
  • 发表时间: 2025年11月11日 (arXiv v2)
  • 论文链接: https://arxiv.org/abs/2404.05969

摘要

本文的主要目标是提供一种随机采样Butcher树的算法，用于常微分方程(ODEs)的概率数值求解。该方法通过消除生成随机分支时间的需求，补充并简化了最近提出的ODE解的概率表示方法。论文通过数值实验比较了树的随机采样与Butcher级数的有限阶截断方法。

研究背景与动机

问题背景

  1. 核心问题: 常微分方程的数值求解是科学计算的基础问题。传统方法使用Butcher级数(基于根树枚举和Taylor展开)来表示ODE解，但高阶树的生成计算代价昂贵。
  2. 重要性:
    • Butcher级数为Runge-Kutta方法提供理论基础
    • 在几何数值积分中有广泛应用
    • 对于复杂非线性ODE，需要更高效的数值方法
  3. 现有方法的局限性:
    • 计算复杂性: 截断Butcher级数需要枚举所有n阶树，计算量随阶数指数增长
    • 固定阶数限制: 传统方法在固定阶数处截断，难以自适应调整精度
    • 先前概率方法的复杂性: 文献20中的方法需要生成随机分支时间序列
  4. 研究动机:
    • 利用Monte Carlo方法通过随机树生成估计Butcher级数
    • 提供精度可随迭代次数提升的替代方案
    • 受Feynman-Kac公式在非线性PDE中应用的启发
    • 简化先前的概率表示方法，移除随机分支时间生成步骤

核心贡献

  1. 直接的随机树生成算法: 提出了基于均匀附着(uniform attachment)的Butcher树随机生成方法，不需要生成随机分支时间，相比文献20更简单直接
  2. 概率表示定理: 建立了ODE解的概率表示公式(定理1): x(t)=E[(t−t0)∣T∣F(T)(x0)(∣T∣∨1)p∣T∣]x(t) = \mathbb{E}\left[\frac{(t-t_0)^{|T|}F(T)(x_0)}{(|T| \vee 1)p_{|T|}}\right] 其中T是随机生成的Butcher树
  3. 半线性ODE扩展: 将方法扩展到半线性ODE x˙(t)=Ax(t)+f(x(t))\dot{x}(t) = Ax(t) + f(x(t))，结合Poisson分布树大小和连续时间Markov链(定理2)
  4. 数值实现与比较: 提供完整的Mathematica代码实现，并通过数值实验验证方法的有效性，比较不同概率分布的性能
  5. 理论分析:
    • 证明了标记树的组合性质(引理3)
    • 推导了最优概率分布(方差最小化)
    • 建立了收敛性条件和矩界

方法详解

任务定义

输入:

  • ODE初值问题: x˙(t)=f(x(t))\dot{x}(t) = f(x(t)), x(t0)=x0∈Rdx(t_0) = x_0 \in \mathbb{R}^d
  • 目标时间点 t>t0t > t_0
  • 光滑函数 f:Rd→Rdf: \mathbb{R}^d \to \mathbb{R}^d

输出:

  • 解 x(t)x(t) 在时间 tt 的近似值

约束条件:

  • ff 的各阶导数有界: ∣∇mf(x0)∣≤C|\nabla^m f(x_0)| \leq C 对所有 m≥0m \geq 0
  • 时间区间限制: t∈[t0,t0+1/C)t \in [t_0, t_0 + 1/C)

Butcher树的基础理论

树的定义与表示

根树 τ=(V,E,∙)\tau = (V, E, \bullet) 由顶点集V、边集E和根节点 ∙\bullet 组成。使用B+操作递归定义:

  • [τ1,…,τm][\tau_1, \ldots, \tau_m] 表示创建新根并连接到 τ1,…,τm\tau_1, \ldots, \tau_m 的根

关键函数:

  1. 基本微分 F:T→C∞(Rd,Rd)F: \mathcal{T} \to C^\infty(\mathbb{R}^d, \mathbb{R}^d):
    • F(∅)=IdF(\emptyset) = \text{Id}, F(∙)=fF(\bullet) = f
    • F(τ)=∇mf(F(τ1),…,F(τm))F(\tau) = \nabla^m f(F(\tau_1), \ldots, F(\tau_m)) 对 τ=[τ1,…,τm]\tau = [\tau_1, \ldots, \tau_m]
  2. 对称性 σ(τ)\sigma(\tau):
    • σ([τ1k1,…,τnkn])=∏i=1nki!σ(τi)ki\sigma([\tau_1^{k_1}, \ldots, \tau_n^{k_n}]) = \prod_{i=1}^n k_i! \sigma(\tau_i)^{k_i}
  3. 树阶乘 τ!\tau!:
    • τ!=∣τ∣∏i=1mτi!\tau! = |\tau| \prod_{i=1}^m \tau_i! 对 τ=[τ1,…,τm]\tau = [\tau_1, \ldots, \tau_m]

Butcher级数表示

ODE解的经典Butcher级数展开: x(t)=∑τ∈T(t−t0)∣τ∣τ!σ(τ)F(τ)(x0)x(t) = \sum_{\tau \in \mathcal{T}} \frac{(t-t_0)^{|\tau|}}{\tau! \sigma(\tau)} F(\tau)(x_0)

系数 α(τ)=∣τ∣!τ!σ(τ)\alpha(\tau) = \frac{|\tau|!}{\tau! \sigma(\tau)} 表示树 τ\tau 的标记方式数量。

标记树与组合结构

标记树定义: 树 τ=(V,E,1)\tau = (V, E, 1) 的顶点用 {1,…,n}\{1, \ldots, n\} 标记，满足父节点标签小于子节点标签。记为 Tn♯\mathcal{T}_n^\sharp。

关键引理(引理3): 任何标记树 τ∈Tn♯\tau \in \mathcal{T}_n^\sharp 可唯一表示为: τ=∙∗l1∙∗l2⋯∗ln−1∙\tau = \bullet *_{l_1} \bullet *_{l_2} \cdots *_{l_{n-1}} \bullet 其中 (l1,…,ln−1)∈△n−1:={(l1,…,ln−1):1≤li≤i}(l_1, \ldots, l_{n-1}) \in \triangle_{n-1} := \{(l_1, \ldots, l_{n-1}): 1 \leq l_i \leq i\}

嫁接积(Grafting product): τ1∗lτ2\tau_1 *_l \tau_2 表示将 τ2\tau_2 的根附着到 τ1\tau_1 的标签为 ll 的顶点。

推论1: nn 阶标记树的数量为 ∣Tn♯∣=(n−1)!|\mathcal{T}_n^\sharp| = (n-1)!

随机树生成算法(算法6)

步骤:

  1. 选择树大小: 从概率分布 (pn)n≥0(p_n)_{n \geq 0} 中采样树的阶数 nn，即 P(∣T∣=n)=pnP(|T| = n) = p_n
  2. 初始化: 从根节点 ∙\bullet (标签1)开始
  3. 迭代附着: 对于 l=1,…,n−1l = 1, \ldots, n-1:
    • 均匀随机选择当前树的一个顶点
    • 将新顶点(标签 l+1l+1)附着到选中的顶点
    • 重复直到达到阶数 nn

条件分布: 给定 ∣T∣=n|T| = n，随机树 TT 在 Tn♯\mathcal{T}_n^\sharp 上均匀分布: qn♯(τ):=P(T=τ∣∣T∣=n)=1(n−1)!q_n^\sharp(\tau) := P(T = \tau | |T| = n) = \frac{1}{(n-1)!}

忽略标记后的条件分布: qn(τ)=P(ι(T)=τ∣∣T∣=n)=α(τ)(n−1)!q_n(\tau) = P(\iota(T) = \tau | |T| = n) = \frac{\alpha(\tau)}{(n-1)!}

概率表示定理

定理1(主要结果): 假设 ∣∇mf(x0)∣≤C|\nabla^m f(x_0)| \leq C 对所有 m≥0m \geq 0，则对 t∈[t0,t0+1/C)t \in [t_0, t_0 + 1/C): x(t)=E[(t−t0)∣T∣F(T)(x0)(∣T∣∨1)p∣T∣]x(t) = \mathbb{E}\left[\frac{(t-t_0)^{|T|}F(T)(x_0)}{(|T| \vee 1)p_{|T|}}\right]

证明思路:

  1. 利用标记树的均匀分布性质
  2. 应用全期望公式: E[⋅]=∑n=0∞pn∑τ∈Tn♯qn♯(τ)F(τ)(x0)\mathbb{E}[\cdot] = \sum_{n=0}^\infty p_n \sum_{\tau \in \mathcal{T}_n^\sharp} q_n^\sharp(\tau) F(\tau)(x_0)
  3. 由 qn♯(τ)=1/(n−1)!q_n^\sharp(\tau) = 1/(n-1)! 和 α(τ)=∣τ∣!/(τ!σ(τ))\alpha(\tau) = |\tau|!/(\tau! \sigma(\tau))，得到Butcher级数
  4. 可积性由矩界控制: E[∣(t−t0)∣T∣F(T)(x0)(∣T∣∨1)p∣T∣∣q]≤∣x0∣qp0q−1+∑n=1∞(C(t−t0))nqnqpnq−1\mathbb{E}\left[\left|\frac{(t-t_0)^{|T|}F(T)(x_0)}{(|T| \vee 1)p_{|T|}}\right|^q\right] \leq \frac{|x_0|^q}{p_0^{q-1}} + \sum_{n=1}^\infty \frac{(C(t-t_0))^{nq}}{n^q p_n^{q-1}}

半线性ODE扩展(定理2)

对于半线性ODE: x˙(t)=Ax(t)+g(x(t))\dot{x}(t) = Ax(t) + g(x(t))，其中 AA 是线性算子:

方法:

  1. 使用Poisson分布生成树大小: pn=e−(t−t0)(t−t0)n/n!p_n = e^{-(t-t_0)}(t-t_0)^n/n!
  2. 引入独立的连续时间Markov链 (Xt)t≥t0(X_t)_{t \geq t_0}，生成元为 AA
  3. 利用指数Butcher级数表示

概率表示: xi(t)=et−t0E[((∣Tt∣−1)∨0)!(Fg(Tt)(x0))Xt−T∣Tt∣1{T∣Tt∣≤t}∣Xt0=i]x_i(t) = e^{t-t_0} \mathbb{E}\left[((|T_t|-1) \vee 0)! (F_g(T_t)(x_0))_{X_{t-T_{|T_t|}}} \mathbf{1}_{\{T_{|T_t|} \leq t\}} \mid X_{t_0} = i\right]

其中 TtT_t 是Poisson大小的随机树，FgF_g 是 gg 的基本微分。

技术创新点

  1. 消除分支时间: 与文献20相比，不需要生成随机时间序列 (Ti)i≥1(T_i)_{i \geq 1}，直接通过均匀附着构造树
  2. 组合等价性: 利用标记树与序列 (l1,…,ln−1)∈△n−1(l_1, \ldots, l_{n-1}) \in \triangle_{n-1} 的双射关系(引理3)，建立了简洁的概率构造
  3. 灵活的分布选择: 框架允许任意概率分布 (pn)n≥0(p_n)_{n \geq 0}，可根据方差优化选择
  4. 半线性结构利用: 通过Markov链处理线性部分，随机树处理非线性部分，实现了结构分解
  5. 理论保证: 提供了明确的收敛条件和矩界，确保Monte Carlo估计的可行性

实验设置

测试方程

例1 (方程27): 指数ODE x˙(t)=ex(t),x(0)=x0\dot{x}(t) = e^{x(t)}, \quad x(0) = x_0 解析解: x(t)=−log⁡(e−x0−t)x(t) = -\log(e^{-x_0} - t)，定义域 t∈[0,e−x0)t \in [0, e^{-x_0})

例2 (方程28): 半线性ODE x˙(t)=tx(t)+x2(t),x(0)=1/2\dot{x}(t) = tx(t) + x^2(t), \quad x(0) = 1/2 解析解: x(t)=et2/22−∫0tes2/2dsx(t) = \frac{e^{t^2/2}}{2 - \int_0^t e^{s^2/2}ds}

实验参数

截断Butcher级数:

  • 阶数范围: n=1,…,8n = 1, \ldots, 8
  • 使用命令 B[f,t,x0,t0,n] 实现

Monte Carlo方法:

  • 几何分布:
    • 参数 p=0.5p = 0.5 或 p=0.75p = 0.75
    • 样本数: 70,000 (方程27), 10,000 (方程28)
  • Poisson分布:
    • 参数 λ=t−t0\lambda = t - t_0
    • 样本数: 100,000
  • 最优分布: p0=c0x0p_0 = c_0 x_0, pn=c0(Ct)n/np_n = c_0(Ct)^n/n (n≥1n \geq 1)

评价指标

  1. 计算时间: 比较不同方法达到相似精度所需的时间
  2. 数值误差: 与解析解的绝对误差
  3. 方差分析: 比较不同概率分布的二阶矩界: x02p0+∑n=1∞(C(t−t0))2nn2pn\frac{x_0^2}{p_0} + \sum_{n=1}^\infty \frac{(C(t-t_0))^{2n}}{n^2 p_n}

实现细节

Mathematica代码:

树生成过程:

  • 使用图结构存储树
  • 顶点标签存储导数信息
  • 随机选择: RandomVariate[DiscreteUniformDistribution[{1, l}]]

实验结果

计算时间比较(表1)

阶数 nn12345678MC (几何)
方程27 (d=1)0s0s0.1s0.1s0.4s0.5s3s21s22s (70K)
方程28 (d=2)0s0s0s0.2s1s13s222s>1h164s (10K)

观察:

  • Butcher级数计算时间随阶数指数增长
  • Monte Carlo方法时间相对稳定
  • 对于方程28，8阶截断超过1小时，而MC方法164秒

主要数值结果(图2)

方程27 (x0=1x_0 = 1, t∈[0,0.35]t \in [0, 0.35]):

  • B-8级数: 与精确解高度吻合
  • B-6级数: 在 t>0.25t > 0.25 处出现偏差
  • MC方法(几何分布, 70K样本): 与精确解吻合良好，方差较小

方程28 (x0=1/2x_0 = 1/2, t∈[0,1]t \in [0, 1]):

  • B-7级数: 精确度高
  • B-5级数: 在 t>0.6t > 0.6 处明显偏离
  • MC方法(几何分布, 10K样本): 跟踪精确解，但方差略大

关键发现:

  1. MC方法在相似计算时间内达到与高阶截断相当的精度
  2. MC方法避免了枚举所有树的组合爆炸问题
  3. 样本数可根据精度需求灵活调整

概率分布比较(图3-4)

二阶矩界分析(图3):

  • 最优分布 pn=c0(Ct)n/np_n = c_0(Ct)^n/n: 在所有 CC 值下给出最小方差界
  • 几何分布 (p=0.5p=0.5): 方差界约为最优分布的2-3倍
  • 几何分布 (p=0.75p=0.75): 方差界更高

数值性能(图4):

  • Poisson分布 (100K样本):
    • 波动显著，方差大
    • 在 t>0.2t > 0.2 处误差增大
    • 理论上方差无界(级数发散)
  • 几何分布 (70K样本):
    • 稳定跟踪精确解
    • 方差有界且较小
    • 在 t∈[0,0.35]t \in [0, 0.35] 上表现优异

结论: 几何分布在实践中表现最佳，平衡了方差和计算效率

树生成示例(图1)

展示了3阶和4阶树的系统生成过程:

  • 3阶树: 2种不同结构
  • 4阶树: 3种主要结构
  • 每个顶点标注对应的导数阶数

相关工作

Butcher级数理论

  1. 经典文献:
    • Butcher (1963, 2016, 2021) 1,2,3: 建立了B-级数代数分析框架
    • Hairer等 (2006) 11: 几何数值积分的标准参考
    • Deuflhard & Bornemann (2002) 10: 科学计算中的ODE方法
  2. 计算实现:
    • Ketcheson & Ranocha (2022) 16: Julia中的B-级数完整实现

概率方法

  1. 分支过程:
    • Skorokhod (1964) 22: 分支扩散过程
    • Vatutin (1993) 23,24: 分支过程与随机树理论
    • Ikeda等 (1968-1969) 15: 分支Markov过程
  2. PDE的概率表示:
    • McKean (1975) 19: Brownian运动在KPP方程中的应用
    • Le Jan & Sznitman (1997) 17: 随机级联与Navier-Stokes方程
    • Dalang等 (2008) 6: 波动方程的Feynman-Kac型公式
    • Henry-Labordère等 (2019) 13: 半线性PDE的分支扩散表示
  3. ODE的概率方法:
    • Penent & Privault (2022) 21: 本文简化的前身工作，需要随机分支时间
    • Nguwi等 (2023) 20: 任意阶导数的完全非线性Feynman-Kac公式

指数积分器

  1. 指数Butcher级数:
    • Hochbruck & Ostermann (2010) 14: 指数积分器综述
    • Luan & Ostermann (2013) 18: 刚性情况的指数B-级数

本文的定位

  • 相比21: 移除随机分支时间，算法更简洁直接
  • 相比20: 专注于ODE而非PDE，提供更高效的实现
  • 相比6,13: 扩展到非线性ODE，使用一般树而非线性链
  • 相比经典方法: 提供Monte Carlo替代方案，避免组合爆炸

结论与讨论

主要结论

  1. 理论贡献:
    • 建立了ODE解的新概率表示(定理1)，基于随机Butcher树
    • 证明了标记树与均匀附着过程的等价性(引理3)
    • 扩展到半线性ODE，结合Poisson过程和Markov链(定理2)
  2. 算法优势:
    • 无需生成随机分支时间，实现更简单
    • 避免高阶树的显式枚举，缓解组合爆炸
    • 精度可通过增加样本数灵活提升
  3. 数值验证:
    • 在测试方程上，MC方法与高阶Butcher级数精度相当
    • 计算时间在高阶情况下显著优于级数截断
    • 几何分布在实践中表现最佳

局限性

  1. 收敛速度:
    • Monte Carlo方法的收敛速度为 O(1/N)O(1/\sqrt{N})，慢于确定性高阶方法
    • 对于低维光滑问题，Runge-Kutta方法仍更高效
    • 论文明确指出: "Monte Carlo估计器无法与经典Runge-Kutta方案竞争"
  2. 适用范围限制:
    • 需要导数有界条件: ∣∇mf(x0)∣≤C|\nabla^m f(x_0)| \leq C
    • 时间区间受限: t∈[t0,t0+1/C)t \in [t_0, t_0 + 1/C)
    • 对于刚性问题或长时间积分，条件可能过于严格
  3. 方差问题:
    • Poisson分布理论方差无界
    • 需要仔细选择概率分布以控制方差
    • 最优分布 pn=c0(Ct)n/np_n = c_0(Ct)^n/n 依赖于未知常数 CC
  4. 高维挑战:
    • 多维ODE的代码实现更复杂(见第7节)
    • 树标记和导数计算的维度依赖性增加
    • 数值实验仅限于1-2维情况
  5. 实验局限:
    • 仅测试了两个简单方程
    • 缺乏与其他概率方法(如20)的直接比较
    • 未探讨自适应采样策略

未来方向

  1. 方法改进:
    • 开发自适应重要性采样策略
    • 研究方差缩减技术(如控制变量)
    • 探索并行化实现
  2. 理论扩展:
    • 放宽导数有界条件
    • 扩展到随机微分方程(SDEs)
    • 研究时间步长自适应策略
  3. 应用领域:
    • 扩展到偏微分方程(PDEs)
    • 应用于高维问题(避免维数灾难)
    • 结合深度学习方法

深度评价

优点

  1. 理论创新性 (★★★★☆):
    • 核心创新: 建立了标记树均匀分布与Butcher级数系数的直接联系，通过引理3的双射关系简化了概率构造
    • 数学严谨: 提供完整的收敛性证明和矩界估计
    • 结构洞察: 半线性ODE的分解(线性部分→Markov链，非线性部分→随机树)体现了深刻的结构理解
  2. 算法简洁性 (★★★★★):
    • 实现简单: 相比文献20,21，算法流程大幅简化
    • 代码可读: Mathematica实现清晰，易于理解和复现
    • 开源共享: 提供GitHub代码库，促进研究可复现性
  3. 数学优雅性 (★★★★★):
    • 嫁接积(grafting product)的引入统一了树操作
    • 概率表示公式(18)形式简洁优美
    • 组合与概率的深度融合
  4. 实验设计 (★★★☆☆):
    • 比较了多种概率分布(Poisson, 几何, 最优)
    • 提供了计算时间和精度的定量分析
    • 方差分析有理论支撑

不足

  1. 实用性有限 (★★☆☆☆):
    • 效率问题: 论文自承"无法与经典Runge-Kutta方案竞争"
    • 适用场景窄: 仅在需要避免树枚举的特殊情况下有优势
    • 参数依赖: 最优分布需要知道常数 CC，实际中难以确定
  2. 实验不足 (★★☆☆☆):
    • 测试案例少: 仅两个简单方程，缺乏复杂系统测试
    • 维度限制: 仅测试1-2维，高维性能未知
    • 缺少对比: 未与其他概率方法(如20)直接比较
    • 误差分析浅: 缺乏详细的误差收敛率分析
  3. 理论局限 (★★★☆☆):
    • 时间区间短: t<t0+1/Ct < t_0 + 1/C 限制了长时间积分
    • 光滑性要求高: 需要所有阶导数有界
    • 方差界粗糙: 矩界(20)可能不够紧
  4. 写作问题 (★★★☆☆):
    • 缺乏对"何时使用此方法"的明确指导
    • 与现有方法的优劣势对比不够充分
    • 某些技术细节(如多维实现)放在附录，影响可读性

影响力评估

  1. 理论贡献 (★★★★☆):
    • 为Butcher级数提供了新的概率视角
    • 连接了组合数学、概率论和数值分析
    • 可能启发其他数值方法的概率化改造
  2. 实用价值 (★★☆☆☆):
    • 当前阶段主要是理论探索
    • 实际应用场景有限
    • 可能在特定问题(如不确定性量化)中有用
  3. 可复现性 (★★★★★):
    • 开源代码完整
    • 算法描述清晰
    • 数值结果可验证
  4. 学术影响:
    • 引用潜力: 中等。方法新颖但实用性限制了应用范围
    • 后续研究: 可能激发方差缩减、自适应采样等改进工作
    • 跨学科价值: 连接了概率、组合、数值分析，有一定跨学科意义

适用场景

推荐使用:

  1. 高阶树枚举困难: 当需要非常高阶Butcher级数且树枚举不可行时
  2. 不确定性量化: 需要同时估计解及其不确定性时
  3. 教学演示: 作为Butcher级数的概率解释工具
  4. 理论研究: 探索数值方法的概率基础

不推荐使用:

  1. 常规ODE求解: Runge-Kutta等经典方法更高效
  2. 实时计算: Monte Carlo方差导致结果不稳定
  3. 刚性问题: 时间步长限制 t<t0+1/Ct < t_0 + 1/C 过于严格
  4. 高精度需求: 收敛速度 O(1/N)O(1/\sqrt{N}) 较慢

综合评分

  • 创新性: 8/10 (新颖的概率视角，简化了先前方法)
  • 严谨性: 9/10 (数学证明完整，理论基础扎实)
  • 实用性: 4/10 (当前阶段实用价值有限)
  • 实验性: 5/10 (实验设计合理但范围有限)
  • 影响力: 6/10 (理论贡献显著，实际应用受限)

总体: 这是一篇理论上优雅、数学上严谨的论文，为Butcher级数提供了全新的概率解释。算法的简洁性和理论的完整性是其主要亮点。然而，实用价值受限于Monte Carlo方法的固有缺陷(收敛慢、方差大)和严格的适用条件。论文更适合作为理论探索和方法论贡献，而非实际求解器的替代方案。未来若能发展出有效的方差缩减技术和自适应策略，该方法的实用性可能显著提升。

参考文献(精选)

  1. Butcher, J.C. (2021). B-Series: Algebraic Analysis of Numerical Methods. Springer. Butcher级数的权威专著
  2. Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric numerical integration. Springer. 几何数值积分经典教材
  3. Penent, G., & Privault, N. (2022). Numerical evaluation of ODE solutions by Monte Carlo enumeration of Butcher series. BIT Numerical Mathematics, 62:1921-1944. 本文简化的前身工作
  4. Henry-Labordère, P., et al. (2019). Branching diffusion representation of semilinear PDEs and Monte Carlo approximation. Ann. Inst. H. Poincaré Probab. Statist., 55(1):184-210. PDE的分支扩散表示
  5. Ketcheson, D.I., & Ranocha, H. (2022). Computing with B-series. ACM Transactions on Mathematical Software. B-级数的Julia实现