Integrate Forced Andronov-Hopf Bifurcation

[iranroman]

Hello torchdiffeq

Thanks for the awesome work.

I am interested in integrating a Hopf oscillator. In fact, I have been able to do this in TensorFlow, but I'm considering a PyTorch port.

Anyways, here's the Hopf oscillator (with hard-coded limit-cycle dynamics after defining values for alpha and beta):

import torch
from torchdiffeq import odeint
import math

class ode_fn(torch.nn.Module):

    def __init__(self, x0, t):
        
        super(ode_fn, self).__init__()        
        
        self.alpha = None
        self.beta = None
        self.x0 = x0
        self.t = t
    
    def forward(self, t, z):
        
        return z * (self.alpha + 1j*2*math.pi + self.beta * torch.pow(torch.abs(z),2))
    
    def solve_ode(self, alpha, beta):
        
        self.alpha = alpha
        self.beta = beta
        
        return odeint(self.forward, self.x0, self.t, method='rk4')

# Test
z0 = torch.tensor(0.01+1j*0).cuda()
t = torch.linspace(0, 10, 1000).cuda()
model = ode_fn(z0, t)

with torch.no_grad():
    z_out = model.solve_ode(alpha = torch.tensor(1).cuda(), beta = -torch.tensor(1).cuda())

Plotting z_out does show a limit-cycle, Yay!

HOWEVER, I would like to integrate a forced Andropov-Hopf oscillator (receiving an input) of the form:

dz/dt = z(alpha + i*2*pi + beta*|z|^2) + x(t)

Note the x(t), which could be another function, OR more specifically in my case, I would like to be able to use time-sampled values in a vector. For example, a time-series with the temperature in Paris yesterday (100 points per hour).

Thank you very much in advance for your help (and sorry if this has been asked before or is somewhere).

[iranroman]

Am I reading correctly that this could potentially be done using something like this example? Or is there a simpler way to do it?

[iranroman]

I came to this solution, but not sure if this is the best way. The oscillator also has a "Hebbian frequency learning" term, which allows the oscillator to match the stimulus frequency. How off am I?

import torch
from torchdiffeq import odeint
import math

def get_force(t,fs,x):
  return x[math.floor(t * fs)]

class ode_fn(torch.nn.Module):

    def __init__(self, init_conds, f0, alpha, beta, cs, cw, cr, x, fs, t):
        
        super().__init__()        
        
        self.alpha = alpha
        self.beta = beta
        self.f0 = f0
        self.cs = cs
        self.cr = cr
        self.cw = cw                
        
        self.init_conds = init_conds
        
        self.t = t        
        self.fs = fs
        self.x = x
    
    def forward(self, t, states):
                
        z,f = states
        w = f * 2 * math.pi

        F = self.cs * get_force(t, self.fs, self.x) * (1 / (1 - torch.conj(z)))        
        freq_learn = (self.cr / torch.abs(z)) * (torch.sin(torch.angle(z))) * F
        elast = self.cw * (w - (self.f0 * 2 * math.pi)) / (self.f0 * 2 * math.pi)    

        z_dot = f * (z * (self.alpha + 1j*2*math.pi + self.beta * torch.pow(torch.abs(z),2)) + F )
        w_dot = f * ( -freq_learn - elast)                

        return [z_dot, w_dot/(2 * math.pi)]
    
    def solve_ode(self):
                
        return odeint(self.forward, self.init_conds, self.t, method='rk4')


# oscillator initial conditions
f0 = 6.73
init_conds = (torch.tensor(0.1+1j*0).cuda(), torch.tensor(f0).cuda())

# oscillator parameters
alpha = torch.tensor(0.233).cuda()
beta = torch.tensor(-2.66).cuda()
cs = torch.tensor(0.51).cuda()
cr = torch.tensor(0.68).cuda()
cw = torch.tensor(0.42).cuda()
f0 = torch.tensor(f0).cuda()

# create forcing signal
fs = 100
dur = 10
t_end = dur - 1/fs
t = torch.linspace(0, 10, dur * fs).cuda()
stim_freq = 5
x = torch.sin(stim_freq * 2 * math.pi * t)

# create model and integrate
model = ode_fn(init_conds, f0, alpha, beta, cs, cw, cr, x, fs, t[:-1])
with torch.no_grad():
    z_out = model.solve_ode()

import matplotlib.pyplot as plt

plt.subplot(3,1,1)
plt.plot(z_out[1].real.cpu().numpy())
plt.subplot(3,1,2)
plt.plot(z_out[0].real.cpu().numpy())
plt.subplot(3,1,3)
plt.plot(x.cpu().numpy())
plt.show()

[iranroman]

I have worked on an implementation that at least behaves as expected. You can check it out here. I would like to know if there is a better way to do this, so feedback welcome! https://github.com/iranroman/DESSEO/blob/main/desseo/models.py