# 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.symbolic import SymbolicSys
sp.init_printing()
%matplotlib inline
print(sp.__version__)
```

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)]
```


```python
y0 = [0, 1]
mu = 2.5
tend = 25
```


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


```python
# Let us plot using 30 data points
res1 = odesys1.integrate(np.linspace(0, tend, 20), y0, [mu], name='vode')
res1.plot()
print(res1.yout.shape)
```


```python
# Let us interpolate between data points
res2 = odesys1.integrate(np.linspace(0, tend, 20), y0, [mu], integrator='cvode', nderiv=1)
res2.plot(m_lim=21)
print(res2.yout.shape)
```


```python
odesys1.integrate(np.linspace(0, tend, 20), y0, [mu], integrator='cvode', nderiv=2)
xplt, yplt = odesys1.plot_result(m_lim=21, interpolate=30)
print(odesys1._internal[1].shape, yplt.shape)
```

Equidistant points are not optimal for plotting this function. Using ``roots`` kwarg we can make the solver report the output where either the function value, its first or second derivative is zero.


```python
odesys2 = SymbolicSys.from_other(odesys1, roots=odesys1.exprs + (odesys1.dep[0],))
# We could also add a higher derivative: tuple(odesys1.get_jac().dot(odesys1.exprs)))
```


```python
# Let us plot using 10 data points
res2 = odesys2.integrate(np.linspace(0, tend, 20), y0, [mu], integrator='cvode',
                                     nderiv=1, atol=1e-4, rtol=1e-4)
xout, yout, info = res2
xplt, yplt = odesys2.plot_result(m_lim=21, interpolate=30, indices=[0])
xroots, yroots = info['roots_output'][0], info['roots_output'][1][:, 0]
plt.plot(xroots, yroots, 'bd')
print(odesys2._internal[1].shape, yplt.shape, xroots.size)
```


```python
odesys2.roots
```


```python
res2.plot(indices=[0])
plt.plot(xplt, [res2.at(_)[0][0, 0] for _ in xplt])
```


```python
res1.plot(indices=[0])
plt.plot(xplt, [res1.at(_, use_deriv=True)[0][0] for _ in xplt])
plt.plot(xplt, [res1.at(_, use_deriv=False)[0][0] for _ in xplt])
```

e6336ccdcfe2988e4d81ef075dbd84dc2ce922c3
ipynb
Jupyter Notebook
examples/van_der_pol_interpolation.ipynb
slayoo/pyodesys
8e1afb195dadf6c6f8e765873bc9dd0fae067c39
BSD-2-Clause
2015-09-29T16:51:03.000Z
2022-02-02T13:26:50.000Z
examples/van_der_pol_interpolation.ipynb
slayoo/pyodesys
8e1afb195dadf6c6f8e765873bc9dd0fae067c39
BSD-2-Clause
2015-09-29T14:40:45.000Z
2021-09-18T19:29:50.000Z
examples/van_der_pol_interpolation.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