Solver for stiff differential equations

New to Julia
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