how do i simulate differential jerk systems in python

Question

I am a beginner in python and I'm trying to simulate a multi scroll chaotic attractor for a project. The details and parameters are in https://sci-hub.se/10.1142/s0218127420501862 [3.2 (i)]

im getting an empty plot. any help would be greatly appreciated. i have tried the odeint for solving the equations.

def system(state, t, alpha):
    x1, x2, x3, x4, x5, x6 = state
    dxdt = [
        x2,
        x3,
        x4,
        x5,
        x6,
        alpha[0] * (f(x1) - x1) - alpha[1] * x2 - alpha[2] * x3 - alpha[3] * x4 - alpha[4] * x5 - alpha[5] * x6
    ]
    return dxdt

# Define the nonlinear function f(x)
def f(x):
    M = 1
    e = 1
    q = 0.01
    return (e / (2 * q)) * sum([np.abs(x - 2 * m * e + q) - np.abs(x - 2 * m * e - q) for m in range(-M, M + 1)])

# Set parameters
alpha_values = [7.5, 10, 10, 40, 10, 10]  
initial_state = [0, 0, 0, 0, 0, 0]
time = np.linspace(0, 100, 10000) 

# Solve ODEs
solution = odeint(system, initial_state, time, args=(alpha_values,))

Answer 1

Well, here you go. The major problem with your code seems to be the initial values (it really doesn't like all zeroes). There is some sensitivity to the parameter R (I found the paper you are alluding to in the end). You also need a solution to a larger Tmax.

Link to paper: https://doi.org/10.1142/S0218127420501862 but I'm relying on university access to that - it may not be an open-access paper.

import numpy as np
from scipy.integrate import odeint

def system(state, t, alpha):
    x1, x2, x3, x4, x5, x6 = state
    dxdt = [
        x2,
        x3,
        x4,
        x5,
        x6,
        alpha[0] * (f(x1) - x1) - alpha[1] * x2 - alpha[2] * x3 - alpha[3] * x4 - alpha[4] * x5 - alpha[5] * x6
    ]
    return dxdt

# Define the nonlinear function f(x)
def f(x):
    M = 1
    e = 1
    q = 0.01
    return (e / (2 * q)) * sum([np.abs(x - 2 * m * e + q) - np.abs(x - 2 * m * e - q) for m in range(-M, M + 1)])

# Set parameters
R = 15.00
Tmax = 1000
Nt = 10000
alpha_values = [7.5, 10.0, R, 40, R, R ]
initial_state = 6 * [0.01]
time = np.linspace(0, Tmax, Nt )

# Solve ODEs
sol = odeint(system, initial_state, time, args=(alpha_values,))

# Plot solution
import matplotlib.pyplot as plt
plt.plot( sol[:, 0], sol[:,1], 'b' )
plt.xlabel( 'x1' )
plt.ylabel( 'x2' )
plt.show()

And here's your multi-scroll attractor. It's quite pretty!