In [9]:
import matplotlib.pyplot as plt
from scipy import *
from scipy import integrate
from scipy.integrate import odeint
import numpy as np
from IPython.display import Image
%matplotlib inline

Lotka Volterra system

$\frac{dx}{dt} = a x -b xy$

$\frac{dy}{dt}= -cy +d xy$

In [10]:
def LotkaVolterra(z,t,a,b,c,d):
    x,y=z
    dydt = [a*x-b*x*y,-c*y+d*x*y]
    return dydt
In [11]:
A=2.0/3.0
B=4.0/3.0
C=1
D=1
y0 = [1.5, 0.5]
t = np.linspace(0, 20, 1001)
sol = odeint(LotkaVolterra, y0, t, args=(A,B,C,D))
In [12]:
fig = plt.figure(num=1,figsize=(10,10))


plt.plot(t, sol[:, 0], 'b', label='x(t) = Prey')
plt.plot(t, sol[:, 1], 'g', label='y(t) = Predator')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
In [20]:
fig = plt.figure(num=1,figsize=(10,10))
ax=fig.add_subplot(111)

xmin=0;xmax=2;ymin=0;ymax=1;
gridx=30;gridy=30;
x=np.array(range(100))/100.;x=xmin+(xmax-xmin)*x;


plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel(r"$x$",fontsize=18)
plt.ylabel(r"$y$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

A=2.0/3.0
B=4.0/3.0
C=1
D=1

#Vector field
X,Y = np.meshgrid( np.linspace(xmin,xmax,gridx),np.linspace(xmin,xmax,gridy) )
dX = A*X-B*X*Y 
dY = -C*Y + D*X*Y

#Normalize arrows
N = np.sqrt(dX**2+dY**2)  
U, V = dX/N, dY/N
ax.quiver( X,Y,U, V,minshaft=5,width=0.005,headwidth=0,color="black")


# Ploting function

#c=-0.5
#f=c*exp(np.arcsin(x))
ax.plot(sol[:,0],sol[:,1])
sol2 = odeint(LotkaVolterra, [1,0.6],t, args=(A,B,C,D))
ax.plot(sol2[:,0],sol2[:,1])

sol3 = odeint(LotkaVolterra, [1,0.95],t, args=(A,B,C,D))
ax.plot(sol3[:,0],sol3[:,1])

plt.show()

/home/ega/.local/lib/python3.5/site-packages/ipykernel/__main__.py:27: RuntimeWarning: invalid value encountered in true_divide

Duffing Equation

$\ddot{x} + A\dot{x}+Bx+Cx^3 = D cos(\omega t)$

Image from Scholarpedia

In [5]:
Image("duffing.png",width=600)
Out[5]:
In [21]:
def LotkaVolterra(z,t,a,b,c,d):
    x,v=z
    dzdt = [v,-a*v-b*x-c*x**3+d*np.cos(t)]
    return dzdt
In [62]:
A=0.2
B=-1
C=1
D=0.3
z0 = [1.5, 0.5]
z1 = [1.5002,0.5]

t = np.linspace(0,200, 1001)
sol1 = odeint(LotkaVolterra, z0, t, args=(A,B,C,D))

fig = plt.figure(num=1,figsize=(10,10))

plt.plot(t, sol1[:, 0], 'b', label='x(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.ylabel('x')
plt.figure(num=1,figsize=(10,10))
plt.show()


t = np.linspace(0, 55, 1001)
sol = odeint(LotkaVolterra, z0, t, args=(A,B,C,D))
sol2 = odeint(LotkaVolterra, z1, t, args=(A,B,C,D))

fig = plt.figure(num=1,figsize=(10,10))

plt.plot(t, sol[:, 0], 'b', label='x(t)')
plt.plot(t, sol2[:, 0], 'r', label='x(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.ylabel('x')
plt.figure(num=1,figsize=(10,10))
plt.show()

fig = plt.figure(num=1,figsize=(10,10))

plt.plot(sol[:, 0],sol[:,1], 'b', label='x(t)')
plt.plot(sol2[:, 0],sol2[:,1], 'r', label='x(t)')

plt.legend(loc='best')
plt.xlabel('x')
plt.ylabel('y')
plt.figure(num=1,figsize=(10,10))
plt.show()

Attractor is a Fractal!

In [7]:
Image("Sierpinski_triangle1.png",width=1000)
Out[7]:
In [ ]: