# Solver for stiff differential equations

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
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)
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
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]


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

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.

this seems not to solve it :-/

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

    dY[1] = sdot
dY[2] = taudot


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
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
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
end


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

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
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
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
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)


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

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

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.

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
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)
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
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)


5 Likes

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
σ = 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)
τ = Y[2]  # shear stress
θ = Y[3]  # state variable
ϵ = 1e6   # small perturbation
ṡ = exp(-μ/a)*v0 * (v0*θ/L)^(-b/a) * exp(τ/a/σ-1) # velocity
τ̇ = k * (vload - ṡ)                               # shear stress
θ̇ = 1 - (ṡ * (θ - ϵ) / L)                         # state evolution
dY[1] = ṡ
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

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