# Solver for stiff differential equations

**URL:** <https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [May 9, 2022, 6:00pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772 "2022-05-09T18:00:10Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![graeffd](https://avatars.discourse-cdn.com/v4/letter/g/ce7236/32.png) [@graeffd](https://discourse.julialang.org/u/graeffd)\
**Post date:** [May 9, 2022, 6:00pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/1 "2022-05-09T18:00:10Z")

</div>

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!

```nohighlight
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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/9/5/95743f389b62d35b5f6fe0223a5869f0fe126db8.png)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 9, 2022, 6:34pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/2 "2022-05-09T18:34:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 9, 2022, 6:36pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/3 "2022-05-09T18:36:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![graeffd](https://avatars.discourse-cdn.com/v4/letter/g/ce7236/32.png) [@graeffd](https://discourse.julialang.org/u/graeffd)\
**Post date:** [May 9, 2022, 6:39pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/4 "2022-05-09T18:39:49Z")

</div>

this seems not to solve it :-/

---

<div class="post-metadata">

**Author:** ![graeffd](https://avatars.discourse-cdn.com/v4/letter/g/ce7236/32.png) [@graeffd](https://discourse.julialang.org/u/graeffd)\
**Post date:** [May 9, 2022, 6:40pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/5 "2022-05-09T18:40:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 9, 2022, 6:41pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/6 "2022-05-09T18:41:41Z")

</div>

```julia
    dY[1] = sdot
    dY[2] = taudot
    dY[3] = thetadot

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 9, 2022, 6:43pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/7 "2022-05-09T18:43:19Z")

</div>

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

```julia
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/](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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 9, 2022, 6:45pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/8 "2022-05-09T18:45:04Z")

</div>

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

```julia
# 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](https://global.discourse-cdn.com/julialang/original/3X/0/d/0d4a27a01f689efb3d62549c3f4fcad247d4d150.png)

---

<div class="post-metadata">

**Author:** ![graeffd](https://avatars.discourse-cdn.com/v4/letter/g/ce7236/32.png) [@graeffd](https://discourse.julialang.org/u/graeffd)\
**Post date:** [May 9, 2022, 6:46pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/9 "2022-05-09T18:46:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 9, 2022, 6:48pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/10 "2022-05-09T18:48:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [May 9, 2022, 7:10pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/11 "2022-05-09T19:10:32Z")

</div>

> [@graeffd](#):
>
> `sb_funcs`

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](https://github.com/JuliaPy/PyCall.jl) to call the Python version of the ODE derivative function and compare it in Julia to what you get from `sb_funcs`.

---

<div class="post-metadata">

**Author:** ![graeffd](https://avatars.discourse-cdn.com/v4/letter/g/ce7236/32.png) [@graeffd](https://discourse.julialang.org/u/graeffd)\
**Post date:** [May 9, 2022, 7:24pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/12 "2022-05-09T19:24:10Z")

</div>

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:

```nohighlight
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](https://global.discourse-cdn.com/julialang/original/3X/5/9/59ec34e19f392854d8a0875e0a726d3663afa435.png)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 10, 2022, 4:21am UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/13 "2022-05-10T04:21:10Z")

</div>

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

```julia
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

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

```

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

---

<div class="post-metadata">

**Author:** ![graeffd](https://avatars.discourse-cdn.com/v4/letter/g/ce7236/32.png) [@graeffd](https://discourse.julialang.org/u/graeffd)\
**Post date:** [May 10, 2022, 3:21pm UTC](https://discourse.julialang.org/t/solver-for-stiff-differential-equations/80772/14 "2022-05-10T15:21:11Z")

</div>

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