Solving a system of ODE using scipy

Question

I am trying to solve the following system of ODE $$ \pmatrix{y_0'\\y_1'} =\pmatrix{-\omega_0&-\omega_1\cos(\omega t)\\-\omega_1\cos(\omega t)&\omega_0} \pmatrix{y_0\\y_1} \tag1$$

This is a system of a spin 1/2 particle subjected to a magnetic field $\omega_0$ and an oscillating field $\omega_1\cos(\omega t)$. The expected result should be a sinusoidal evolution of $y_0$ and $y_1$.

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

w0 = 1e3
w1 = 5e3
w=1e3
args = [w0,w1,w]
y0, t0 = [1,0],0 #initial values
def func(t,y,args):
    w0,w1,w = args

    f=[-w0*y[0]-w1*np.cos(w*t)*y[1],
       -w1*np.cos(w*t)*y[0]+w0*y[1]]
    return f

r = ode(func).set_integrator('zvode', method='bdf')
r.set_initial_value(y0,t0).set_f_params(args)
t = np.linspace(0,100e-3,100)
dt = t[1]-t[0]
y=[]
while r.successful() and r.t < t[-1]:
    y.append(r.integrate(r.t+dt))
plt.plot(np.array(y)[:,0])
plt.plot(np.array(y)[:,1])

Using scipy's generic ODE solver, the results seem to blow up ($y_0\approx e^{140}$). I am not sure If I understood how to translate (1) into the function func.

Answer 1

The eigenvalues of the system matrix are on average $\pm\sqrt{ω_0^2+\frac12ω_1^2}\sim 3.5\cdot 10^3$. The roughly estimated size of the solution at time $t=0.1$ is $\exp(3.5\cdot 10^3\cdot 0.1)=e^{350}\sim 10^{152}$. This is qualitatively conform with what you observed.

As a spin system you want the spinor to stay on the unit ball in $\Bbb C^2$ which is achieved theoretically if the system matrix in $y'=Ay$ is skew-symmetric, $A^*=-A$. To turn the symmetric matrix skew-symmetric, you need to multiply it by $i$.