# Van der Pol oscillator
We will look at the second order differentual equation (see https://en.wikipedia.org/wiki/Van_der_Pol_oscillator):

$$
{d^2y_0 \over dx^2}-\mu(1-y_0^2){dy_0 \over dx}+y_0= 0
$$


```python
from __future__ import division, print_function
import itertools
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
#from pyodesys.native.gsl import NativeGSLSys as SymbolicSys
from pyodesys.native.cvode import NativeCvodeSys as SymbolicSys
sp.init_printing()
%matplotlib inline
print(sp.__version__)
```

Note that we imported ``NativeCvodeSys`` as ``SymbolicSys``, this speed up the time of integration by more than an order of magnitude due to using compiled C++ code for our mathematical expressions.

One way to reduce the order of our second order differential equation is to formulate a system of first order ODEs, using:

$$ y_1 = \dot y_0 $$

which gives us:

$$
\begin{cases}
\dot y_0 = y_1 \\
\dot y_1 = \mu(1-y_0^2) y_1-y_0
\end{cases}
$$

Let's call this system of ordinary differential equations vdp1:


```python
vdp1 = lambda x, y, p: [y[1], -y[0] + p[0]*y[1]*(1 - y[0]**2)]
mu_val = 2.5
y0_1 = [0.0, 1.0]
y0_1, (y0_1[0], vdp1(0, y0_1, [mu_val])[0])
```

An alternative would be to use use the Liénard transformation:

$$ y = x - x^3/3 - \dot x/\mu $$


```python
transf = lambda y, dydx, p: [y[0], y[0] - y[0]**3/3 - dydx[0]/p[0]]
x, mu = sp.symbols('x mu', real=True)
y = [yi(x) for yi in sp.symbols('y:2', cls=sp.Function)]
dydx = [yi.diff(x) for yi in y]
[sp.Eq(yi, expr, evaluate=False) for yi, expr in zip(y, transf(y, dydx, [mu]))]  # Just for displaying
```

which gives us (we could generate this result using SymPy):
$$
\begin{cases}
\dot y_0 = \mu \left(y_0-\frac{1}{3}y_0^3-y_1\right) \\
\dot y_1 = \frac{1}{\mu} y_0
\end{cases}
$$


```python
vdp2 = lambda x, y, p: [p[0]*(y[0] - y[0]**3/3 - y[1]), y[0]/p[0]]
calc_y0_2 = lambda y0, mu: transf(y0, vdp1(0, y0, [mu]), [mu])
y0_2 = calc_y0_2(y0_1, mu_val)
(y0_2, y0_2[0], vdp2(0, y0_2, [mu_val])[0])
```


```python
def solve_and_plot(odesys, y0, tout, mu, indices=None, integrator='native', **kwargs):
    plt.figure(figsize=(16, 4))
    xout, yout, info = odesys.integrate(tout, y0, [mu], integrator=integrator, **kwargs)
    plt.subplot(1, 2, 1)
    odesys.plot_result(indices=indices, ls=('-',), c=('k', 'r'))
    plt.legend(loc='best')
    plt.subplot(1, 2, 2)
    odesys.plot_phase_plane()
    info.pop('internal_xout')  # too much output
    info.pop('internal_yout')
    return len(xout), info
```


```python
tend = 25
```


```python
odesys1 = SymbolicSys.from_callback(vdp1, 2, 1, names='y0 y1'.split())
odesys1.exprs
```


```python
for mu in [0, 3, 9]:
    solve_and_plot(odesys1, y0_1, np.linspace(0, tend, 500), mu)
```

As we see, the period ($\tau$) varies with $\mu$, in 1952 Mary Cartwright derived an approximate formula for $\tau$ (valid for large $\mu$):

$$
\tau  = (3 - 2 \ln 2)\mu + 2 \alpha \mu^{-1/3}
$$

where $\alpha \approx 2.338$ 


```python
tau = lambda mu: 1.6137056388801094*mu + 4.676*mu**(-1./3)
for mu in [20, 40, 60]:
    solve_and_plot(odesys1, y0_1, np.linspace(0, 5*tau(mu), 500), mu)
```

For larger values of $\mu$ we run into trouble (the numerical solver fails).
The phase portrait is not well resolved due to rapid variations in y1.

Let us look at our alternative formulation:


```python
odesys2 = SymbolicSys.from_callback(vdp2, 2, 1, names='y0 y1'.split())
odesys2.exprs
```


```python
solve_and_plot(odesys2, y0_2, tend, mu_val, nsteps=2000)
```

This looks much better. Let's see if the solver has an easier time dealing with this formulation of y2 for large values of $\mu$:


```python
ls = itertools.cycle(('-', '--', ':'))
for mu in [84, 160, 320]:
    y0_2 = calc_y0_2(y0_1, mu)
    print(y0_2)
    solve_and_plot(odesys2, y0_2, np.linspace(0, 5*tau(mu), 500), mu)
```

Indeed it has.

## Stiffness

Let us compare the performance of explicit and implicit steppers (Adams and BDF respectivecly) for varying values of $\mu$


```python
solve_and_plot(odesys2, calc_y0_2(y0_1, mu_val), tend, mu_val, nsteps=2000)
J = odesys2.get_jac()
J
```

For this simple system we can afford calculating the eigenvalues analytically


```python
odesys2._NativeCode._written_files
```


```python
symbs = odesys2.dep + tuple(odesys2.params)
symbs
```


```python
Jeig = J.eigenvals().keys()
eig_cbs = [sp.lambdify(symbs, eig, modules='numpy') for eig in Jeig]
Jeig
```


```python
eigvals = np.array([(eig_cbs[0](*(tuple(yvals)+(mu_val + 0j,))),
                     eig_cbs[1](*(tuple(yvals)+(mu_val + 0j,)))) for yvals in odesys2._internal[1]])
```


```python
plt.plot(odesys2._internal[0], odesys2.stiffness(), label='from SVD')
plt.plot(odesys2._internal[0], np.abs(eigvals[:,0])/np.abs(eigvals[:,1]), label='analytic')
plt.legend()
```

## Audio
Plotting is instructive from a mathematical standpoint, but these equations were often investigated by to *listening* to audio amplified by electrical circuits modelled by the equation. So let us generate some audio.


```python
def arr_to_wav(arr, rate=44100):
    from IPython.display import Audio
    from scipy.io.wavfile import write
    scaled = np.int16(arr/np.max(np.abs(arr)) * 32767)
    write('test.wav', rate, scaled)
    return Audio('test.wav')
```


```python
xout, yout, info = odesys2.integrate(np.linspace(0, 500*tau(40.0), 2*44100), y0_1, [40.0], integrator='native')
arr_to_wav(yout[:, 0])
```


```python
def overlay(tend_mu, odesys=odesys2, time=3, rate=44100, plot=False):
    yout_tot = None
    for tend, mu in tend_mu:
        xout, yout, info = odesys.integrate(np.linspace(0, tend*tau(mu[0]), time*rate), y0_1, mu, integrator='native')
        print(tend, mu, tend*tau(mu[0]))
        if yout_tot is None:
            yout_tot = yout[:, 0]
        else:
            yout_tot += yout[:, 0]
    if plot:
        plt.figure(figsize=(16,4))
        plt.plot(yout_tot[slice(None) if plot is True else slice(0, plot)])
    return arr_to_wav(yout_tot, rate=rate)
```


```python
overlay([
    (400, [2.0]),
    (410, [2.1]),
], plot=10000)
```

## Forced van der pol oscillator


```python
vdp_forced = lambda x, y, p: [y[1], p[1]*sp.sin(p[2]*x) - y[0] + p[0]*y[1]*(1 - y[0]**2)]
odesys_forced = SymbolicSys.from_callback(vdp_forced, 2, 3)
overlay([(700, [8, 1, 0.5])], odesys_forced, plot=5000)  # Non-chaotic behavior
```


```python
overlay([(700, [8, 1.2, 0.6])], odesys_forced, plot=5000)  # Chaotic behavior
```

### Transient $\mu$


```python
vdp_transient = lambda x, y, p: [y[1], - y[0] + p[0]*sp.exp(-p[1]*x)*y[1]*(1 - y[0]**2)]
odesys_transient = SymbolicSys.from_callback(vdp_transient, 2, 2)
odesys_transient.exprs
```


```python
overlay([
    (440, [0.1, 1/2500.]),
    (445, [0.5, 1/1000.]),
    (890, [0.1, 2/2500.]),
    (896, [0.5, 2/1000.]),
], odesys_transient, plot=-1)
```


```python
odesys2._native._written_files
```


```python

```

7b7635ad2042f3be99512a7785aba26dc9269728
ipynb
Jupyter Notebook
examples/van_der_pol_formulations.ipynb
slayoo/pyodesys
8e1afb195dadf6c6f8e765873bc9dd0fae067c39
BSD-2-Clause
2015-09-29T16:51:03.000Z
2022-02-02T13:26:50.000Z
examples/van_der_pol_formulations.ipynb
slayoo/pyodesys
8e1afb195dadf6c6f8e765873bc9dd0fae067c39
BSD-2-Clause
2015-09-29T14:40:45.000Z
2021-09-18T19:29:50.000Z
examples/van_der_pol_formulations.ipynb
slayoo/pyodesys
8e1afb195dadf6c6f8e765873bc9dd0fae067c39
BSD-2-Clause
2016-03-18T14:00:39.000Z
2021-09-17T13:54:29.000Z
Qwen/Qwen-72B
1. YES
2. YES