import numpy as np;
import matplotlib.pyplot as plt;
plt.rcParams.update({"text.usetex":True});
%config InlineBackend.figure_format = "svg"
from ipywidgets import interactive
from scipy.integrate import solve_ivp
from mayavi import mlab
def mu_effect(mu = 1):
x = np.linspace(-3, 3, 100); y = np.linspace(-3, 3, 100);
X, Y = np.meshgrid(x, y);
U = mu-X**2;
V = -Y;
#plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
plt.streamplot(X, Y, U, V, integration_direction='forward', density=1.5)
ax = plt.gca(); ax.set_aspect(1);
plt.axvline(0); plt.axhline(0);
plt.plot(mu**0.5, 0, 'ok');
plt.plot(-mu**0.5, 0, 'sr');
w = interactive(mu_effect, mu=(-1, 1, 0.1))
w
def snb(a = 0.4, b = 1.0):
x = np.linspace(-3, 3, 100); y = np.linspace(-3, 3, 100);
X, Y = np.meshgrid(x, y);
U = -a*X+Y;
V = -X**2/(1+X**2)-b*Y;
#plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
plt.streamplot(X, Y, U, V, integration_direction='forward', density=1.5)
ax = plt.gca(); ax.set_aspect(1);
plt.axvline(0); plt.axhline(0);
plt.contour(X, Y, U, levels=[0.0])
plt.contour(X, Y, V, levels=[0.0])
#plt.xlim(-2, 0); plt.ylim(-2, 0)
w = interactive(snb, a = (-2, 2, 0.1), b = (-2, 2, 0.1));
w
def mu_effect(mu = 1):
x = np.linspace(-3, 3, 100); y = np.linspace(-3, 3, 100);
X, Y = np.meshgrid(x, y);
U = mu*X-Y+X*Y**2;
V = X+mu*Y+Y**3;
#plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
plt.streamplot(X, Y, U, V, integration_direction='forward', density=1.5)
ax = plt.gca(); ax.set_aspect(1);
plt.axvline(0); plt.axhline(0);
w = interactive(mu_effect, mu=(-1, 1, 0.1))
w
plt.rcParams["figure.figsize"] = [4,4]
def mysys(t, x, gamma, d, omega): # returns the RHS
return ([x[1], x[0]-x[0]**3-gamma*x[1]+d*np.cos(omega*t)]);
gamma = 0.1; d = 0.35; omega = 1.4; tp = 2*np.pi/omega;
Nmax = 30000;
tmax = Nmax*tp
tspan = [0, tmax]
x0 = 0;
y0 = 0;
ics = [x0, y0]
sol = solve_ivp(mysys, tspan, ics, args=(gamma, d, omega), dense_output=True, atol = 1e-8, rtol=1e-8);
tout = np.linspace(0, np.max(tspan), 10000);
xout = sol.sol(tout)[0];
yout = sol.sol(tout)[1];
plt.rcParams["figure.figsize"] = [6,12]
fig, ax = plt.subplots(2, 1)
tlast = np.linspace((Nmax-50)*tp, tmax, 2000);
xlast = sol.sol(tlast)[0];
ylast = sol.sol(tlast)[1];
ax[0].plot(xlast, ylast)
ts = tp*np.arange(Nmax-20000, Nmax)
xs = sol.sol(ts)[0];
ys = sol.sol(ts)[1];
ax[1].plot(xs, ys, '.k', markersize=1)
[<matplotlib.lines.Line2D at 0x1ef175fb310>]
plt.rcParams["figure.figsize"] = [4,4]
def mysys(t, x, s, r, b): # returns the RHS
return ([s*(x[1]-x[0]), r*x[0]-x[1]-x[0]*x[2], x[0]*x[1]-b*x[2]]);
tmax = 200;
tspan = [0, tmax];
s = 10; b = 8/3; r = 28;
x0 = 0;
y0 = 1;
z0 = 0;
ics = [x0, y0, z0]
sol = solve_ivp(mysys, tspan, ics, args=(s, r, b), dense_output=True, atol = 1e-8, rtol=1e-8);
tout = np.linspace(0, np.max(tspan), 20000);
xout = sol.sol(tout)[0];
yout = sol.sol(tout)[1];
zout = sol.sol(tout)[2];
plt.plot(xout, zout);
r = (xout**2+yout**2+zout**2)**0.5
mlab.plot3d(xout, yout, zout, r, tube_radius = 0.2);
mlab.show();
plt.rcParams["figure.figsize"] = [7,14]
from scipy.signal import find_peaks
peaks,_ = find_peaks(zout)
fig, ax = plt.subplots(2,1);
ax[0].plot(tout, zout)
ax[0].plot(tout[peaks], zout[peaks], 'bx')
z = zout[peaks];
zk = z[0:-2]; zkp1 = z[1:-1];
ax[1].plot(zk, zkp1, '.k', markersize=2)
ax[1].set_aspect(1);
xd = np.linspace(30, 45);
ax[1].plot(xd, xd, '-m')
[<matplotlib.lines.Line2D at 0x1ef5f66af70>]
plt.rcParams["figure.figsize"] = [4,4]
def mysys(t, x, s, r, b): # returns the RHS
return ([s*(x[1]-x[0]), r*x[0]-x[1]-x[0]*x[2], x[0]*x[1]-b*x[2]]);
def ycross(t, x, s, r, b):
return x[1]
tmax = 200;
tspan = [0, tmax];
s = 10; b = 8/3; r = 28;
x0 = 0;
y0 = 1;
z0 = 0;
ics = [x0, y0, z0]
sol = solve_ivp(mysys, tspan, ics, events=ycross, args=(s, r, b), dense_output=True, atol = 1e-8, rtol=1e-8);
tout = np.linspace(0, np.max(tspan), 20000);
xout = sol.sol(tout)[0];
yout = sol.sol(tout)[1];
zout = sol.sol(tout)[2];
plt.plot(xout, zout);
te = sol.t_events[0];
xe = sol.y_events[0][:,0]
ye = sol.y_events[0][:,1]
ze = sol.y_events[0][:,2]
r = (xout**2+yout**2+zout**2)**0.5
mlab.plot3d(xout, yout, zout, (yout>0).astype(float), tube_radius = 0.2);
mlab.points3d(xe, ye, ze, scale_factor=0.8)
mlab.show();
