我正在尝试实现 Rackauckas & Nie (2017) 的名为 ESRK1 的 SDE 求解器和名为 RSwM1 的自适应步长算法。我正在写一个python实现，主要是为了向自己确认我已经正确理解了算法。然而，我在实现 ESRK1 时已经遇到了一个问题：当我在描述几何布朗运动的简单 SDE 上以越来越短的时间步长测试我的实现时，随着
dt我相信这个算法是作为 Julia 库 DifferentialEquations.jl 的一部分实现的，所以我想也许我可以通过查看 Julia 代码找到一些帮助。但是，我在查找相关代码时遇到了一些麻烦。如果有人可以向我指出这些算法的相关 Julia librar(y/ies) 中的 ESRK1 和 RSwM1 实现（或者实际上任何其他可读且正确的实现），我将不胜感激。
我在 StochasticDiffEq.jl 的 github 存储库中搜索了 ESRK 和 RSwM，但没有找到任何我可以真正识别为我正在阅读的论文中的方法的内容：
https://github.com/search?q=repo%3ASciML%2FStochasticDiffEq.jl+rswm&type=code
为了完整起见，这是我自己在 python 中尚未正确的 ESRK1 实现：
def ESRK1(U, t, dt, f, g, dW, dZ):
# Implementation of ESRK1, following Rackauckas & Nie (2017)
# Eq. (2), (3) and (4) and Table 1
# Stochastic integrals, taken from Eqs. (25) - (30) in Rackauckas & Nie (2017)
I1 = dW
I11 = (I1**2 - dt) / 2
I111 = (I1**3 - 3*dt*I1) / 6
I10 = (I1 + dZ/np.sqrt(3))*dt / 2
# Coefficients, taken from Table 1 in Rackauckas & Nie (2017)
# All coefficients not included below are zero
c0_2 = 3/4
c1_2, c1_3, c1_4 = 1/4, 1, 1/4
A0_21 = 3/4
B0_21 = 3/2
A1_21 = 1/4
A1_31 = 1
A1_43 = 1/4
B1_21 = 1/2
B1_31 = -1
B1_41, B1_42, B1_43 = -5, 3, 1/2
alpha1, alpha2 = 1/2, 2/3
alpha_tilde1, alpha_tilde2 = 1/2, 1/2
beta1_1, beta1_2, beta1_3 = -1, 4/3, 2/3
beta2_1, beta2_2, beta2_3 = -1, 4/3, -1/3
beta3_1, beta3_2, beta3_3 = 2, -4/3, -2/3
beta4_1, beta4_2, beta4_3, beta4_4 = -2, 5/3, -2/3, 1
# Stages in the Runge-Kutta approximation
# Eqs. (3) and (4) and Table 1 in Rackauckas & Nie (2017)
# First stages
H0_1 = U # H^(0)_1
H1_1 = U
# Second stages
H0_2 = U + A0_21 * f(t, H0_1)*dt + B0_21 * g(t, H1_1)*I10/dt
H1_2 = U + A1_21 * f(t, H0_1)*dt + B1_21 * g(t, H1_1)*np.sqrt(dt)
# Third stages
H0_3 = U
H1_3 = U + A1_31 * f(t, H0_1) * dt + B1_31 * g(t, H1_1) * np.sqrt(dt)
# Fourth stages
H0_4 = U
H1_4 = U + A1_43 * f(t, H0_3) * dt + (B1_41 * g(t, H1_1) + B1_42 * g(t+c1_2*dt, H1_2) + B1_43 * g(t+c1_3*dt, H1_3)) * np.sqrt(dt)
# Construct next position
# Eq. (2) and Table 1 in Rackauckas & Nie (2017)
U_ = U + (alpha1*f(t, H0_1) + alpha2*f(t+c0_2*dt, H0_2))*dt \
+ (beta1_1*I1 + beta2_1*I11/np.sqrt(dt) + beta3_1*I10/dt ) * g(t, H1_1) \
+ (beta1_2*I1 + beta2_2*I11/np.sqrt(dt) + beta3_2*I10/dt ) * g(t + c1_2*dt, H1_2) \
+ (beta1_3*I1 + beta2_3*I11/np.sqrt(dt) + beta3_3*I10/dt ) * g(t + c1_3*dt, H1_3) \
+ (beta4_4*I111/dt ) * g(t + c1_4*dt, H1_4)
# Calculate error estimate
# Eq. (9) and Table 1 in Rackauckas & Nie (2017)
E = -dt*(f(t, H0_1) + f(t + c0_2*dt, H0_2))/6 \
+ (beta3_1*I10/dt + beta4_1*I111/dt)*g(t, H1_1) \
+ (beta3_2*I10/dt + beta4_2*I111/dt)*g(t + c1_2*dt, H1_2) \
+ (beta3_3*I10/dt + beta4_3*I111/dt)*g(t + c1_3*dt, H1_3) \
+ (beta4_4*I111/dt)*g(t + c1_4*dt, H1_4)
# Return next position and error
return U_, E
Rackauckas 和 Nie（2017）：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5844583/pdf/nihms920388.pdf
所以，我想我找到了问题的答案：我试图实现的方法的代码（或至少是 ESRK1 部分）似乎可以在 StochasticDiffEq.jl 的两个地方找到：
论文中的系数（尽管它们现在被称为 SRIW1 而不是 ESRK1）：
以及方法（也可以与其他系数一起使用）：
它的可读性不是很好，至少对于像我这样的非 Julia 程序员来说不是，因为它使用了一些高级功能，但我想我可以稍微耐心地解决它。