In [1]:

    

%run Steppers.ipynb #Import stepper functions (Euler, RK4, etc)
%run ODEs.ipynb     #Import ODEs (Pendulum etc)
import matplotlib
from mpl_toolkits.mplot3d import axes3d
matplotlib.use('nbagg')
import matplotlib.pyplot as plt


    

In [2]:

    

x=np.linspace(-3*np.pi, 3*np.pi,1000)
plt.plot(x,NonlinearSin(x,0))
plt.grid()
plt.show()


    

In [3]:

    

t0 = 0
tf = 7
h  = 0.001
n=int((tf-t0)/h)
t  = np.linspace(t0,tf,n)  #Generate timesteps

for x0 in np.linspace(-6,6,8):
    x = RK4(NonlinearSin, (x0,), t)
    plt.plot(t,x)
plt.grid()
plt.show()


    

In [4]:

    

for r in (1,-1):
    x=np.linspace(0, 4.5, 1000)
    for k in range(1,4):
        plt.plot(x,LogisticEquation(x,0,r=r,k=k), label="k={}".format(k))
    plt.title("Logistic equation r={}".format(r))
    plt.grid()
    plt.legend()
    plt.xlabel("x")
    plt.ylabel(r"$\frac{dx}{dt}$")
    plt.show()


    

In [5]:

    

t0 = 0
tf = 100
h  = 0.005
t  = np.arange(t0,tf,h)  #Generate timesteps

x0=11
r=0.4
k=1

f0=(x0,)
args=(r,k)
for k in [2,4,6,8,10,13]:
    args=(r,k)
    x=RK4(Budworm, f0, t, args)
    plt.plot(t,x,label="k={}".format(k))
plt.grid()
plt.legend()
plt.title("Budworm population dynamics r = {}".format(r))
plt.show()


    

In [6]:

    

def LHS(x,r,k):
    return r*(1-x/k)
def RHS(x,r,k):
    return (x)/(1+x**2)

r=0.4
k=10

x=np.linspace(0,15,1000)

lhs=LHS(x,r,k)
rhs=RHS(x,r,k)

plt.plot(x,lhs, label=r"$r(1-\frac{x}{k})$")
plt.plot(x,rhs, label=r"$\frac{x}{1+x^2}$")
plt.ylim(0)
plt.xlim(0)
plt.xlabel("x")
plt.legend()
plt.show()

x=np.linspace(0,1.1*k,1000)
lhs=LHS(x,r,k)
rhs=RHS(x,r,k)

total = x*(lhs-rhs)

plt.plot(x,total)
plt.xlabel("x")
plt.ylabel(r"$\frac{dx}{dt}$")
plt.grid()
plt.show()

args=(r,k)
for x0 in np.linspace(0,6,7):
    x=RK4(Budworm, (x0,), t, args)
    plt.plot(t,x,label=r"$x_0={:0.2f}$".format(x0))
plt.grid()
plt.title("Budworm population dynamics r = {}, k = {}".format(r,k))
plt.legend()
plt.show()


    

In [36]:

    

def IntersectionIndices(f,g):
    return np.where(np.diff(np.sign(f-g)))[0]

rMin, rMax = (0,1)
kMin, kMax = (0.1,20)
resolution=500

rs=np.linspace(rMin,rMax,resolution)
ks=np.linspace(kMin,kMax,resolution)
points=np.zeros((len(rs),len(ks)))
for i, k in enumerate(ks):
    x=np.linspace(0,k,500)
    for j, r in enumerate(rs):
        lhs=LHS(x,r,k)
        rhs=RHS(x,r,k)
        points[j,i]=len(IntersectionIndices(lhs,rhs))     

plt.contourf(ks, rs, points)
plt.colorbar()
plt.title("Bistable region")
plt.xlabel("k")
plt.ylabel("r")
plt.show()


    

In [ ]: