I'm trying to implement the SDE solver called ESRK1 and the adaptive stepsize algorithm called RSwM1 from Rackauckas & Nie (2017). I'm writing a python implementation, mainly to confirm to myself that I've understood the algorithm correctly. However, I'm running into a problem already at the implementation of ESRK1: When I test my implementation with shorter and shorter timesteps on a simple SDE describing geometric Brownian motion, the solution does not converge as dt
becomes smaller, indicating that I have a mistake in my code.
I believe this algorithm is implemented as part of the library DifferentialEquations.jl in Julia, so I thought perhaps I could find some help by looking at the Julia code. However, I have had some trouble locating the relevant code. If someone could point me to the implementation of ESRK1 and RSwM1 in the relevant Julia librar(y/ies) (or indeed any other readable and correct implementation) of these algorithms, I would be most grateful.
I searched for ESRK and RSwM in the github repo of StochasticDiffEq.jl, but I didn't find anything I could really recognise as the method from the paper I'm reading:
https://github.com/search?q=repo%3ASciML%2FStochasticDiffEq.jl+rswm&type=code
Update: I found the code for ESRK1, as shown in my answer below, but I'm still unable to find the code for RSwM1.
For completeness, here is my own not-yet-correct implementation of ESRK1 in python:
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
source https://stackoverflow.com/questions/77508717/where-to-find-the-code-for-esrk1-and-rswm1-in-the-julia-library-source-code
Comments
Post a Comment