[Python] Solving a set of coupled ordinary differential equation representing the coupled Chua...

I am trying to write a python code to solve the ODE representing a coupled Chua oscillators. The set of ODEs can be found here. The non-linear term is defined as f(x) And the possible solution is given here.

For two Chua oscillators coupled bidirectionally via the (x)-variable, the ODEs are modified to include coupling terms. Let ((x_1, y_1, z_1)) represent the state variables of the first oscillator and ((x_2, y_2, z_2)) represent the state variables of the second oscillator. The coupled system can be written as:

For the first Chua oscillator:


\frac{dx_1}{dt} = \alpha \left( y_1 - x_1 - f(x_1) \right) + \kappa (x_2 - x_1)
\frac{dy_1}{dt} = x_1 - y_1 + z_1
\frac{dz_1}{dt} = -\beta y_1

For the second Chua oscillator:


\frac{dx_2}{dt} = \alpha \left( y_2 - x_2 - f(x_2) \right) + \kappa (x_1 - x_2)
\frac{dy_2}{dt} = x_2 - y_2 + z_2
\frac{dz_2}{dt} = -\beta y_2


Here:

• (\kappa) is the coupling strength between the two oscillators,
• (f(x)) is the same nonlinear function for both oscillators, typically defined as:

f(x) = m_1 x + 0.5 (m_0 - m_1) \left( |x + 1| - |x - 1| \right)


where (m_0) and (m_1) are parameters of the piecewise-linear Chua diode characteristic.

The coupling term (\kappa (x_2 - x_1)) introduces bidirectional coupling between the oscillators.

The article is available on: https://www.researchgate.net/profile/Makoto_Itoh/publication/263913998_Chaos_synchronization_in_Chua's_Circuit/links/58a7e92192851cf0e3b97b0f/Chaos-synchronization-in-Chuas-Circuit.pdf?_tp=eyJjb250ZXh0Ijp7ImZpcnN0UGFnZSI6InB1YmxpY2F0aW9uIiwicGFnZSI6InB1YmxpY2F0aW9uIn19

I have written the code as follows:

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

# Define the nonlinear function f for Chua's circuit
def f(u, m0, m1):
return m1*u + (m0 - m1)/2 * (np.abs(u + 1) - np.abs(u - 1))

# Define the ODE system for two coupled Chua oscillators
def Chua_dr(t, X, ep, alpha, beta, m0, m1):
x, y, z, xp, yp, zp = X
w = ep * (x - xp) # Coupling term
return [
alpha*(y - x - f(x, m0, m1)) - w, # First oscillator equations
x - y + z,
-beta * y,
alpha*(yp - xp - f(xp, m0, m1)) + w, # Second oscillator equations
xp - yp + zp,
-beta * yp
]

# Solve the ODE system
def sol_dr(t, T, ep, alpha, y0_dr):
sol = solve_ivp(
Chua_dr, [0, T], y0_dr,
args=(ep, alpha, beta, m0, m1), method='LSODA',
dense_output=True, rtol=1e-9, atol=1e-12)
return sol.sol(t)

# Parameters for Chua oscillators
alpha, beta, m0, m1 = 10, 14.87, -1.27, -0.68
ep = 0.15 # Coupling strength

# Initial conditions for the system
y0_dr = [-2, 0.2, 4, 0.7, 0.4, -0.8]
# Generate time array
T = 500
t = np.linspace(100, T, (T-100)*100+1)

# Solve the system and plot the results
sol = sol_dr(t, T, ep, alpha, y0_dr)


The numerical solution diverges for coupling value ep>0.12. I have tried few things, such as changing integratin method for solve_ivp from default RK45 to LSODA suitable for stiff ODEs, integrating for shorter time span and using the last value as the initial condition to further the integration ahead. But all gives the same result. The concerned article on the other hand shows results with nice attractor with synchronised manifold for all values of ep.

Continue reading...