Skip to main content

Where to find the code for ESRK1 and RSwM1 in the Julia library source code?

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

Popular posts from this blog

Confusion between commands.Bot and discord.Client | Which one should I use?

Whenever you look at YouTube tutorials or code from this website there is a real variation. Some developers use client = discord.Client(intents=intents) while the others use bot = commands.Bot(command_prefix="something", intents=intents) . Now I know slightly about the difference but I get errors from different places from my code when I use either of them and its confusing. Especially since there has a few changes over the years in discord.py it is hard to find the real difference. I tried sticking to discord.Client then I found that there are more features in commands.Bot . Then I found errors when using commands.Bot . An example of this is: When I try to use commands.Bot client = commands.Bot(command_prefix=">",intents=intents) async def load(): for filename in os.listdir("./Cogs"): if filename.endswith(".py"): client.load_extension(f"Cogs.{filename[:-3]}") The above doesnt giveany response from my Cogs ...

How to show number of registered users in Laravel based on usertype?

i'm trying to display data from the database in the admin dashboard i used this: <?php use Illuminate\Support\Facades\DB; $users = DB::table('users')->count(); echo $users; ?> and i have successfully get the correct data from the database but what if i want to display a specific data for example in this user table there is "usertype" that specify if the user is normal user or admin i want to user the same code above but to display a specific usertype i tried this: <?php use Illuminate\Support\Facades\DB; $users = DB::table('users')->count()->WHERE usertype =admin; echo $users; ?> but it didn't work, what am i doing wrong? source https://stackoverflow.com/questions/68199726/how-to-show-number-of-registered-users-in-laravel-based-on-usertype

Where and how is this Laravel kernel constructor called? [closed]

Where and how is this Laravel kernel constructor called? public fucntion __construct(Application $app, $Router $roouter) { } I have read the documentation and some online tutorial but I can find any clear explanation. I am learning Laravel and I am wondering where does this kernel constructor receives its arguments from. "POSTMOTERM" CLARIFICATION: Here is more clarity.I have checked the boostrap/app.php and it is only used for boostrapping the interfaces into the container class. What is not clear to me is where and how the Kernel class is instatiated and the arguments passed to the object calling the constructor.Something similar to; obj = new kernel(arg1,arg2) or, is the framework using some magic functions somewhere? Special gratitude to those who burn their eyeballs and brain cells on this trivia before it goes into a full blown menopause alias "MARKED AS DUPLICATE". To some of the itchy-finger keyboard warriors, a.k.a The mods,because I believe in th...