I was investigating the following system of two weakly coupled identical Van der Pol oscillators
$$\left\{\begin{array}{@{}l@{}} \ddot{x}_1 + x_1 + \epsilon(x_1^2 - 1)\dot{x}_1 = \epsilon k(x_2 - x_1)\\ \\ \ddot{x}_2 + x_2 + \epsilon(x_2^2 - 1)\dot{x}_2 = \epsilon k(x_1 - x_2) \end{array}\right.\,$$
As I understand it, for small epsilon, I can treat both oscillators as perturbed harmonic oscillators. So supposing that $x_i$ and $\dot{x}_i$ take the form: $$x_i = r_i(t)\cos(t+\phi_i(t)) \qquad \dot{x}_i = r_i(t)\sin(t+\phi_i(t))$$ Or in exponential form, with $z_i(t) = \frac{r_i}{2}e^{i\phi_i}$: $$x_i = z_i(t)e^{it} + z_i(t)^* e^{-it} \qquad \dot{x}_i = i(z_i(t)e^{it}-z_i(t)^*e^{-it})$$
with $r_i(t)$ and $\phi_i(t)$ being slowly varying functions of time.
Substituting these into the coupled system. $$2i\dot{z}_2(t)e^{it} - z_2(t)e^{it} - z_2(t)^* e^{-it} + \left(z_2(t)e^{it} + z_2(t)^* e^{-it} \right) \\ + i\epsilon \left[ z_2^3(t) e^{i3t} - z_2^2(t)z(t)^* e^{it} + z_2^2(t)^* z_2(t)e^{-it} - z_2^3(t)^* e^{-i3t} + 2z_2^2(t)z(t)^* e^{it} - 2z_2(t)z_2^2(t)^* e^{-it} - z_2(t)e^{it} + z_2(t)^* e^{-it} \right] \\ = \epsilon k\left( z_1(t)e^{it} + z_2(t)^*e^{-it} - z_2(t)e^{it} - z_2(t)^*e^{-it} \right)$$
Multiplying by $e^{-it}$ then performing first order averaging, I obtain the following slow flow :
$$ \left\{\begin{array}{@{}l@{}} \dot{r}_1 = \frac{\epsilon k}{2} r_2 \sin(\phi_2 - \phi_1) + \frac{\epsilon}{8} r_1 \left( 4 - r_1^2 \right)\\ \\ \dot{r}_2 = \frac{\epsilon k}{2} r_1 \sin(\phi_1 - \phi_2) + \frac{\epsilon}{8} r_2 \left( 4 - r_2^2 \right)\\ \\ \dot{\phi}_1 = \frac{\epsilon k}{2}\left( 1 - \frac{r_2}{r_1}\cos(\phi_2 - \phi_1)\right) \\ \\ \dot{\phi}_2 = \frac{\epsilon k}{2}\left( 1 - \frac{r_1}{r_2}\cos(\phi_1 - \phi_2)\right) \end{array}\right.\,$$ I was interested in comparing how faithful of an approximation this was, so I plugged both the slow flow system and the original system of two oscillators into a numerical integrator. Using the following python code to compute and graph the solutions of the original and the averaged system.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.markers import MarkerStyle
from scipy.integrate import solve_ivp
# y = dx/dt
# random initial conditions between -0.5 and 0.5, [x1_0, y1_0, x2_0, y2_0]
x0 = ((np.random.rand(1, 4)[0])-0.5)
print(x0)
# Converting to initial conditions to polar coordinates for slow flow
z1 = x0[0] + x0[1]*1j
print(np.abs(z1), np.rad2deg(np.angle(z1)))
z2 = x0[2] + x0[3]*1j
r0 = [np.abs(z1), np.abs(z2), np.angle(z2)-np.angle(z1)]
print(r0)
# [r1_0, r2_0, phi1_0, phi2_0]
r_phi0 = [np.abs(z1), np.abs(z2), np.angle(z1), np.angle(z2)]
print("compare initial conditions")
print(f"[{r_phi0[0]*np.cos(r_phi0[2])}, {r_phi0[0]*np.sin(r_phi0[2])}, {r_phi0[1]*np.cos(r_phi0[3])}, {r_phi0[1]*np.sin(r_phi0[3])}]")
EPSILON = 0.05
KAPPA = 1
t_start = 0
t_stop = 500
DURATION = t_stop
SAMPLE_RATE = 1000
SAMPLES = int(DURATION*SAMPLE_RATE)
t = np.linspace(t_start, t_stop, SAMPLE_RATE*DURATION)
# numerical integration of original system
def coupled_vdp_deriv(t, X, eps, k):
x1, y1, x2, y2 = X
return [y1, -eps*(x1*x1 - 1)*y1 - x1 + eps*k*(x2 - x1), y2, -eps*(x2*x2 - 1)*y2 - x2 + eps*k*(x1 - x2)]
sol = solve_ivp(coupled_vdp_deriv,
[t_start, t_stop],
x0,
args=[EPSILON, KAPPA],
method="LSODA",
dense_output=True,
t_eval=t)
x1_ori, y1_ori, x2_ori, y2_ori = sol.y
# numerical integration of averaged system
def coupled_r_phi(t, R, eps, k):
r1, r2, phi1, phi2 = R
return [ eps/8*r1*(4-r1**2) + eps*k*r2/2*np.sin(phi2-phi1), eps/8*r2*(4-r2**2) - eps*k*r1/2*np.sin(phi2-phi1), eps*k/2*(1 - r2/r1*np.cos(phi2 - phi1)), eps*k/2*(1 - r1/r2*np.cos(phi2 - phi1)) ]
sol = solve_ivp(coupled_r_phi,
[t_start, t_stop],
r_phi0,
args=[EPSILON, KAPPA],
method="LSODA",
dense_output=True,
t_eval=t)
r1, r2, phi1, phi2 = sol.y
# reconstructing the 'fast flow' variable
x1_av = r1*np.cos(t + phi1)
x2_av = r2*np.cos(t + phi2)
# visualising the data side by side
fig, axes = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(20, 8))
#axes[0].set_ylim(-0.3, 10)
axes[0].set_ylabel("$x_original$")
axes[0].set_xlabel("$t$")
axes[0].plot(t, x1_ori, label="$x_1$", color="blue", alpha=0.8)
axes[0].plot(t, x2_ori, label="$x_2$", color="red", alpha=0.8)
axes[0].legend()
axes[1].set_ylabel("$x_averaged$")
axes[1].plot(t, x1_av, label="$x_1$", color="blue", alpha=0.8)
axes[1].plot(t, x2_av, label="$x_2$", color="red", alpha=0.8)
axes[1].legend()
dirpath="figs/coupled_comparison"
fig.suptitle(f"$X_0$=[{x0[0]:.3f}, {x0[1]:.3f}, {x0[2]:.3f}, {x0[3]:.3f}] $k$={KAPPA} $\epsilon$={EPSILON}")
plt.show()
I was expecting to find discrepancies between the two such as small errors in amplitude or phase. However, I found instead that the two often displayed completely different behaviour from the get go (wildly different amplitudes), while converging towards the same limit behaviour.
So I'm wondering what this could be? I checked to make sure that the initial conditions given to the integrator were equivalent by converting the polar ones back to cartesian coordinates. And they were (except for tiny floating point errors at the 9th decimal digit).
I know the forced Van der Pol system can exhibit chaos with certain driving frequencies. And I know floating point errors in chaotic systems can grow tremendously, but the behaviour here is different, even at the very beginning.
Or is the averaged system simply not faithful enough to the original system? And if so, why?
| Example 1 | Example 2 | Example 3 |
|---|---|---|
 |
 |
 |
Solution
Nothing weird happening. As @Lutz Lehmann pointed out,there was an orientation problem due to a missing negative sign in the code, compounded by a transcription error in the post for the trigonometric form of $\dot{x}$: $$\dot{x}_i = -r_i(t)\sin(t+\phi_i(t))$$ The initial conditions in polar form weren't correct:
# Converting to initial conditions to polar coordinates for slow flow
z1 = x0[0] - x0[1]*1j
z2 = x0[2] - x0[3]*1j
You have a common problem of orientation that is built-in into the harmonic equation $\ddot x+x=0$. The usual state parametrization is as $u=(u_1,u_2)=(x,\dot x)$. The derivative is then $\dot u=(\dot x,\ddot x)=(u_2,-u_1)$.
Now however you want to switch to polar coordinates $u=A(\cos\theta,\sin\theta)$ with the same component order, so that $x=A\cos\theta$, $\dot x=A\sin \theta$. However, the derivative is now $A\dot\theta(-\sin\theta,\cos\theta)$. As one can see, the orientation of the rotation is opposite to the one wanted above (with the intuitive assumption of $\dot\theta>0$; the actual solution has $\dot\theta=-1$, which reverses the orientation).
One easy fix is just to accept that $\theta=\theta_0-t$. Another variant is to set $u_2=-\dot x$, so that $\dot u_1=-u_2$, $\dot u_2=u_1$ reflects the derivative structure of the polar coordinates. As far as I can see, the second variant is the minimal change for your code, at the start you have to set $z_i=x_i-\mathfrak i \,\dot x_i$ to get the polar coordinates right, everything else remaining unchanged. With that the plots you use become visually almost the same.