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 $$
In [ ]:
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:
In [ ]:
vdp1 = lambda x, y, p: [y[1], -y[0] + p[0]*y[1]*(1 - y[0]**2)]
In [ ]:
y0 = [0, 1]
mu = 2.5
tend = 25
In [ ]:
odesys1 = SymbolicSys.from_callback(vdp1, 2, 1, names='y0 y1'.split())
odesys1.exprs
In [ ]:
# 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)
In [ ]:
# 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)
In [ ]:
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.
In [ ]:
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)))
In [ ]:
# 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)
In [ ]:
odesys2.roots
In [ ]:
res2.plot(indices=[0])
plt.plot(xplt, [res2.at(_)[0][0, 0] for _ in xplt])
In [ ]:
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])