Example forward model of rate-and-state friction in python

rsfs.py

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
from math import exp,log

def friction(t,y):
    k = 1e-3
    a = 0.005
    b = 0.01
    Dc = 10.
    mu_ref = 0.6
    V_ref = 1.

    if t < 10:
        V_lp = 1.
    else:
        V_lp = 10.

    # Just to help readability
    #y[0] is mu (friction)
    #y[1] is theta

    n = len(y)
    dydt = np.zeros((n,1))

    v = V_ref * exp((y[0] - mu_ref - b * log(V_ref * y[1] / Dc)) / a)
    dydt[0] = k * (V_lp - v)
    dydt[1] = 1. - v * y[1] / Dc
    print v
    return dydt

r = integrate.ode(friction).set_integrator('vode', order=5,max_step=0.001,method='bdf',atol=1e-10,rtol=1e-6)

k = 1e-3
a = 0.005
b = 0.01
Dc = 10.
mu_ref = 0.6

# Time range
t_start = 0.0
t_final = 50.
delta_t = 0.01
num_steps = np.floor((t_final-t_start)/delta_t)+1

# Initial conditions
mu_t_zero = 0.6
V_ref = 1.
theta_t_zero = Dc/V_ref
v = V_ref
r.set_initial_value([mu_t_zero, theta_t_zero], t_start)

# Create arrays to store trajectory
t = np.zeros((num_steps,1))
mu = np.zeros((num_steps,1))
theta = np.zeros((num_steps,1))
velocity = np.zeros((num_steps,1))
t[0] = t_start
mu[0] = mu_ref
theta[0] = theta_t_zero
velocity[0] = v

# Integrate the ODE(s) across each delta_t timestep
k = 1
while r.successful() and k < num_steps:
    #integrate.ode.set_f_params(r,velocity,k)
    r.integrate(r.t + delta_t)

    # Store the results to plot later
    t[k] = r.t
    mu[k] = r.y[0]
    theta[k] = r.y[1]
    k += 1

# Make some plots

cjm_t,cjm_mu = np.loadtxt('cjm_step.tim',skiprows=2,unpack=True)

fig = plt.figure()
ax1 = plt.subplot(311)
ax1.plot(t, mu,color='r')
ax1.set_xlim(t_start, t_final)
ax1.plot(cjm_t+10.,cjm_mu,color='k')
ax1.set_xlabel('Time [sec]')
ax1.set_ylabel('Friction')
ax1.grid('on')

ax2 = plt.subplot(312)
ax2.plot(t, theta, 'r')
ax2.set_xlim(t_start, t_final)
ax2.set_xlabel('Time [sec]')
ax2.set_ylabel('State Variable')
ax2.grid('on')

ax3 = plt.subplot(313)
ax3.plot(t, velocity, 'r')
ax3.set_xlim(t_start, t_final)
ax3.set_xlabel('Time [sec]')
ax3.set_ylabel('Velocity')
ax3.grid('on')

plt.show()