Solver for stiff differential equations

New to Julia
1

Hi everybody,

I’m now to Julia and trying to translate a simple slider-block model from Python. The code is modelling earthquake cycles. It consists mainly of a set of 3 differential equations. In Python, I use from scipy.intergrate the function solve_ivp which solves the differential equations with many different solvers reliably. I figured out that the solver ‘lsoda’ is fasted for what I need.
However, on my translated code in Julia, where I use ‘solve’ also with ‘lsoda’, the equations are not solved correctly. Essentially the solver misses the earthquake cycles, which are very short. I expect a behaviour like in the attached plot. However, Julia does not replicate, what I would expect.

Here is a minimal code example. I’m happy for any help!

using DifferentialEquations
using LSODA

#set parameters here
a = 1.899e-2        # material property
b = 2.959e-2        # material property
mu = 5e-02          # steady-state friction (no steady state at mu=0)
sigma = 1.5e5       # Pa normal stress
k = 1e8             # N/m spring constant
vload = 1.6e-7      # m/s starting velocity
v0 = 2.2e-6         # m/s reference velocity 
L = 9.0e-6 ;        # m characterisitc slip distance

# set of differential equations
function sb_funcs(dY,Y,pars,t)
    a,b,mu,sigma,k,vload,v0,L = pars
    tau = Y[2]   # shear stress
    theta = Y[3] # state variable
    eps = 1e6    # small perturbation
    sdot = exp(-mu/a)*v0 * (v0*theta/L)^(-b/a) * exp(tau/a/sigma-1) # velocity
    taudot = k * (vload - sdot)                                     # shear stress
    thetadot = 1 - (sdot * (theta - eps) / L)                       # state evolution
    dY = [sdot, taudot, thetadot]
    return dY
end

#run the model
tspan = (0.0,7200)                     # time span
pars_in = (a,b,mu,sigma,k,vload,v0,L)  # input parameters
y0 = [0., 1e-6, 1e-6]                  # starting values
#define the ODE Problem
prob = ODEProblem(sb_funcs,y0,tspan,pars_in)
#solve differential equations
sol = solve(prob, lsoda(), abstol=1e-9, reltol=1e-9, save_everystep=true)

The solver returns:

retcode: Success
Interpolation: 1st order linear
t: 7-element Vector{Float64}:
    0.0
    0.22768399153212332
    0.45536798306424664
 2277.2952833042973
 4554.13519862553
 6830.9751139467635
 7200.0
u: 7-element Vector{Vector{Float64}}:
 [0.0, 1.0e-6, 1.0e-6]
 [0.0, 1.0e-6, 1.0e-6]
 [-1.1281505e-317, 1.0e-6, 1.0e-6]
 [-1.2854371873e-314, 1.0e-6, 1.0e-6]
 [1.2845366769083356e-306, 1.0e-6, 1.0e-6]
 [2.927253055270818e-303, 1.0e-6, 1.0e-6]
 [2.1337310042281397e-302, 1.0e-6, 1.0e-6]

example_sliderblock
example_sliderblock1800×1200 94.6 KB

2

You could perhaps try setting dtmax=0.1 (or some appropriate value) when calling solve to see if the solver catches the events properly.

3

Actually, the problem is likely that you fail to write intibdY and instead create a new dY. Try writing the derivative into the array that is sent into the dynamics function.

4

this seems not to solve it :-/

5

I do not understand entirely, what you mean. Could you give an example?

6 
    dY[1] = sdot
    dY[2] = taudot
    dY[3] = thetadot
7

For better performance, I’d move the constants into the dynamics function or declare them as constants, i.e.

function sb_funcs(dY,Y,pars,t)
    #set parameters here
    a = 1.899e-2        # material property
    b = 2.959e-2        # material property
    mu = 5e-02          # steady-state friction (no steady state at mu=0)
    sigma = 1.5e5       # Pa normal stress
    k = 1e8             # N/m spring constant
    vload = 1.6e-7      # m/s starting velocity
    v0 = 2.2e-6         # m/s reference velocity 
    L = 9.0e-6 ;        # m characterisitc slip distance
    a,b,mu,sigma,k,vload,v0,L = pars
    tau = Y[2]   # shear stress
    theta = Y[3] # state variable
    eps = 1e6    # small perturbation
    sdot = exp(-mu/a)*v0 * (v0*theta/L)^(-b/a) * exp(tau/a/sigma-1) # velocity
    taudot = k * (vload - sdot)                                     # shear stress
    thetadot = 1 - (sdot * (theta - eps) / L)                       # state evolution
    dY[1] = sdot
    dY[2] = taudot
    dY[3] = thetadot
end

for additional performance tips, see
https://docs.julialang.org/en/v1/manual/performance-tips/

EDIT: constants were passed as parameters and the performance penalty from globals does thus not apply

8

Are you sure the dynamics implemented in julia is the same as in python?
I get from

# set of differential equations
function sb_funcs(dY,Y,pars,t)
    #set parameters here
    a = 1.899e-2        # material property
    b = 2.959e-2        # material property
    mu = 5e-02          # steady-state friction (no steady state at mu=0)
    sigma = 1.5e5       # Pa normal stress
    k = 1e8             # N/m spring constant
    vload = 1.6e-7      # m/s starting velocity
    v0 = 2.2e-6         # m/s reference velocity 
    L = 9.0e-6 ;        # m characterisitc slip distance
    a,b,mu,sigma,k,vload,v0,L = pars
    tau = Y[2]   # shear stress
    theta = Y[3] # state variable
    eps = 1e6    # small perturbation
    sdot = exp(-mu/a)*v0 * (v0*theta/L)^(-b/a) * exp(tau/a/sigma-1) # velocity
    taudot = k * (vload - sdot)                                     # shear stress
    thetadot = 1 - (sdot * (theta - eps) / L)                       # state evolution
    dY[1] = sdot
    dY[2] = taudot
    dY[3] = thetadot
end

#run the model
tspan = (0.0,7200)                     # time span
pars_in = (a,b,mu,sigma,k,vload,v0,L)  # input parameters
y0 = [0., 1e-6, 1e-6]                  # starting values
#define the ODE Problem
prob = ODEProblem(sb_funcs,y0,tspan,pars_in)
#solve differential equations
sol = solve(prob, lsoda(), abstol=1e-9, reltol=1e-9, save_everystep=true)
plot(sol, layout=3)

image
image859×629 37.2 KB

9

cool! it is starting to do something. The output does not seem correct, yet, but I think this is already one step closer!

10

BTW, sorry for the mistaken suggestion about the global constants, I see that you pass them as parameters, and you will thus not suffer performance penalty from initially having them as globals.

1 Like
11

Don’t debug via the ODE solvers, just directly call sb_funcs in Julia and Python and confirm you are getting the same result to floating point accuracy. You can use PyCall.jl to call the Python version of the ODE derivative function and compare it in Julia to what you get from sb_funcs.

12

Oh, no problem!
I also found two typos in the code. Two parenthesis and a minus sign got lost in the code translation from python to Julia. NOW IT WORS!
Thank you so much! I hope that implementing this into my more complex versions of the model will speed up my code by a factor of 10-100 compared to Python.

For completeness, here is the correctly working code:

using DifferentialEquations
using LSODA
using Plots

#set parameters here
a = 1.899e-2        # material property
b = 2.959e-2        # material property
mu = 5e-02          # steady-state friction (no steady state at mu=0)
sigma = 1.5e5       # Pa normal stress
k = 1e8             # N/m spring constant
vload = 1.6e-7      # m/s starting velocity
v0 = 2.2e-6         # m/s reference velocity 
L = 9.0e-6 ;        # m characterisitc slip distance

# set of differential equations
function sb_funcs(dY,Y,pars,t)
    a,b,mu,sigma,k,vload,v0,L = pars
    tau = Y[2]   # shear stress
    theta = Y[3] # state variable
    eps = 1e-6   # small perturbation
    sdot = exp(-mu/a)*v0 * (v0*theta/L)^(-b/a) * (exp(tau/a/sigma)-1) # velocity
    taudot = k * (vload - sdot)                                     # shear stress
    thetadot = 1 - (sdot * (theta - eps) / L)                       # state evolution
    dY[1] = sdot
    dY[2] = taudot
    dY[3] = thetadot
    return dY
end

#run the model
tspan = (0.0,7200)                     # time span
pars_in = (a,b,mu,sigma,k,vload,v0,L)  # input parameters
y0 = [0., 1e-9, 1e-9]                  # starting values
#define the ODE Problem
prob = ODEProblem(sb_funcs,y0,tspan,pars_in)
#solve differential equations
sol = solve(prob, lsoda(), abstol=1e-9, reltol=1e-9)
plot(sol, layout=3)

Screenshot 2022-05-09 at 12.23.46

5 Likes
13

Since you’re new to julia and use a lot of greek letters typed out, I just figure’d I’d show you how pretty the code could look with unicode variables

a = 1.899e-2        # material property
b = 2.959e-2        # material property
μ = 5e-02           # steady-state friction (no steady state at μ=0)
σ = 1.5e5           # Pa normal stress
k = 1e8             # N/m spring constant
vload = 1.6e-7      # m/s starting velocity
v0 = 2.2e-6         # m/s reference velocity 
L = 9.0e-6 ;        # m characterisitc slip distance

# set of differential equations
function sb_funcs(dY,Y,pars,t)
    a,b,μ,σ,k,vload,v0,L = pars
    τ = Y[2]  # shear stress
    θ = Y[3]  # state variable
    ϵ = 1e6   # small perturbation
    ṡ = exp(-μ/a)*v0 * (v0*θ/L)^(-b/a) * exp(τ/a/σ-1) # velocity
    τ̇ = k * (vload - ṡ)                               # shear stress
    θ̇ = 1 - (ṡ * (θ - ϵ) / L)                         # state evolution
    dY[1] = ṡ
    dY[2] = τ̇
    dY[3] = θ̇
end

whether or not you prefer this is of course up to you, some people do, some do not. If you ever see a unicode symbol in code and want to know how to type it, simply paste it into the REPL in help mode

help?> θ
"θ" can be typed by \theta<tab>

(Note: the code above is the old version without your fixes of typos)

4 Likes
14

Looks indeed much easier to read with the greek symbols and the time derivatives represented by the dots! Thanks for showing this to me!

1 Like