0
$\begingroup$

I am recently working on the code based on the stick-slip phenomenon in Python. It's the stick-slip oscillator (chapter 3.4) in https://www.sciencedirect.com/science/article/pii/S0888327020301205. And my aim is to achieve the plot as in the Fig. 9. The problem is that after some calculations I keep getting the following error:

Traceback (most recent call last):
  File "stick_slip.py", line 89, in <module>
    cur_y = list(sol.y_events[-1][0])
IndexError: index 0 is out of bounds for axis 0 with size 0

And I have no idea with what the problem is, so I would be very grateful if you could help me with that error and eliminate it from my code or tell what should I correct here to get expected results.

The code is below. So I've started with inclusion of libraries, then definition of parameters, initial conditions, etc.

#libraries, etc.
import scipy as sp
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint 
from scipy.integrate import solve_ivp
#define variables and parameters
eta=2.088
v0=1
fs=8
fk=5.0
epsilon = 1e-12

period = (2.0*np.pi)/eta
transient_periods = 40
periods_to_save = 100
points_per_period = 100
tau_transient=[*np.linspace(0, transient_periods*period, transient_periods*points_per_period+1)]
tau_to_save = [*np.linspace(transient_periods*period, (transient_periods+periods_to_save)*period, points_per_period*periods_to_save+1)]
y=[0,0,0]
# t_span = (0,period,100000)
#tau_2=(period,2000,1000)

Then I set stick phenomenon as a global variable.

is_stick = False

Next, I defined function based on the conditions in article. So what are the faces of the equations when stick exist, what if not

def f(tau, y):
    dy3 = eta
    if is_stick:
        dy1=v0
        dy2=0
    elif np.abs(y[2] - v0) == 0:
        fff = fk*np.sign(y[0]-np.cos(y[3])) 
        dy1 = y[2]
        dy2 = -y[0]+np.cos(y[3])+fff
    else:
        fff = fk*np.sign(v0-y[2])
        dy1 = y[2]
        dy2 = -y[0]+np.cos(y[3])+fff
    
    return [dy1, dy2, dy3]     

Then I used scipy.integrate.solve_ivp to make some events, so to stop integrating one function when the event occurs, and then moving to the other function and so on.

#from slip to stick, detect slip and move to stick   
def sliptostick(tau, y): return y[2]-v0 
sliptostick.terminal=True 
sliptostick.direction=0
   

# detect a stick and move to slip 
def sticktoslip(tau, y): return abs(y[0] - np.cos(y[3])) - fs
sticktoslip.terminal=True
sticktoslip.direction=1
   

Here I set some necessary values.

cur_tau = 0.0
cur_y = y[:]
tau_max = tau_transient[-1]

tau_list = []
y_list = []

Then is the part that update the stick and value of y.

def update_stick(y):
    global is_stick
    if np.abs(v0-y[2]) == 0.0 and abs(y[0]-np.cos(y[3]))-fs < 0:
        is_stick = True
    else:
        is_stick = False

update_stick(cur_y)

print(cur_tau)
print(tau_max)
#input()

Now it's the part of the code with the while loop, so here, the program is integrating stick or slip funcion (depending on which event ocurred). The data are saving to appropriate variables.

while tau_transient:
    
    if(is_stick):
        sol = solve_ivp(f,(cur_tau, tau_max),y,events=sticktoslip, t_eval = tau_transient)
    else:
        print("cur_tau: " + str(cur_tau))
        print("cur_y: " + str(cur_y))
        print("tau_max: " + str(tau_max))
        print("tau_trans_len:" + str(len(tau_transient)))
        #input()
        sol = solve_ivp(f,(cur_tau, tau_max),y,events=sliptostick,t_eval = tau_transient)

    tau_list.extend([*sol.t])
    y_list.extend([[*q] for q in list((sol.y).T)])
    cur_y = list(sol.y_events[-1][0])
    cur_tau = sol.t_events[-1][0]  
    update_stick(cur_y)
    tau_transient = [local_t for local_t in tau_transient if local_t > cur_tau]

Here, I plot the results.


plt.plot(tau_list, [q[0] for q in y_list],'k-',label='Y')
plt.plot(tau_list, [q[2] for q in y_list],'r-',label='dY/dt')
plt.xlabel('Time')
plt.legend()
plt.grid()
plt.show()

plt.plot([q[0] for q in y_list],[q[2] for q in y_list],'dodgerblue',markersize=1)
plt.xlabel('y1')
plt.ylabel('y2')
plt.grid()
plt.show()

$\endgroup$
3
  • 1
    $\begingroup$ I believe this question is more suited for StackOverflow. Still, it is likely that the integration reaches the end time tau_max without having triggered any events, hence the list of occurrences of your event is empty. Try and look at the trajectory you obtain to see if that makes sense. Then maybe you need to continue the integration until you have an event that is triggered. $\endgroup$
    – Laurent90
    Mar 29 at 20:59
  • 1
    $\begingroup$ Print sol.message and evaluate sol.status, 0=finished, 1=terminal event. // It should not have any influence here, but in general sol.t_events[-1][0] addresses the first occurrence of the last event in the event list. $\endgroup$ Mar 29 at 21:45
  • $\begingroup$ I think that adding the ODE to your question might help potential answers. $\endgroup$
    – nicoguaro
    Mar 30 at 18:55
0
$\begingroup$

I believe this does what you are looking for:

import math
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

#--------------------------------------
# Problem parameters
#--------------------------------------
# --- coefs ---------------------------
eta = 2.088
v0 = 1.
fs = 8.
fk = 5.
# --- time ----------------------------
period = (2.*np.pi) / eta
n_periods = 400
t_final = n_periods*period
max_timestep = 1e-3
# --- initial values ------------------
t = 0.
x = np.zeros(3)

#--------------------------------------
# ODE functions
#--------------------------------------
def f(x):
    """small f function."""
    if not math.isclose(x[1], v0):
        return fk * np.sign(v0 - x[1])

    value = x[0] - np.cos(x[2])

    if np.abs(value) >= fs:
        return fk * np.sign(value)

    return value


def fun(t, x):
    """ODE function dx/dt (t) = fun(t, x(t))."""
    return np.array([x[1], -x[0] + np.cos(x[2]) + f(x), eta])

#--------------------------------------
# Event detection
#--------------------------------------
def stickevent(t, x): 
    """detect sticking event."""
    return x[1]-v0

stickevent.terminal = True 


def slipevent(t, x): 
    """detect slip event."""
    return np.abs(x[0] - np.cos(x[2])) - fs

slipevent.terminal = True

def stick_phase(x):
    """check if we are sticking."""
    if math.isclose(x[1], v0):
        return True
    return False

#--------------------------------------
# Main loop
#--------------------------------------
t_grid = []
x_vals = []

max_iter = 1000
iter_count = 0

while t < t_final and iter_count < max_iter:
    
    if stick_phase(x): event = slipevent
    else: event = stickevent
    print(f"It: {iter_count+1} ### Sticking: {stick_phase(x)} ### Time: {t} of {t_final}", end='\r')
    
    sol = solve_ivp(fun, (t, t_final), x, events=event, max_step=max_timestep)

    t_grid.extend(sol.t)
    x_vals.extend((np.squeeze(y) for y in np.hsplit(sol.y, sol.y.shape[1])))
    
    t = t_grid[-1]
    x = x_vals[-1]

    iter_count += 1


#--------------------------------------
# Plotting
#--------------------------------------
x1 = [x[0] for x in x_vals]
x2 = [x[1] for x in x_vals]

plt.plot(x1, x2, color='k')
plt.xlabel('x1')
plt.ylabel('x2')
plt.savefig('stick-slip.jpg')

It produces:

Plot

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.