Problem 1

<br data-mce-bogus="1">
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

#numerical parameters
n = 31 #length of array
dt = 10. #time step

#Verhulst parameters
alpha = .1 #growth rate
K = 1000. #carrying capacity

#create arrays
N_numerical = np.zeros((n,3)) #make N arrray
N_analyical = np.zeros((n,1)) #make N arrray
t = np.linspace(0,(float(n)-1.)*dt,num=n) #make t array at constant interval

#initial condition
N_initial = 1.0

###FUNCTIONS
#Verhulst Analytical Solution
def verhulst_analytical(alpha,K,N_initial,t):
    N = (K * N_initial * np.exp(alpha*t))/(K+N_initial*(np.exp(alpha*t)-1.))
    return N

#Verhulst Numerical Integrator
def verhulst_rate(N,alpha,K):
    rate = alpha * N * (1.0 - N / K)
    return rate

#Euler
def euler(Nini,tini,alpha,K,dt,n):
    N = np.zeros((n))
    N[0] = Nini
    for i in xrange(1,n):
        rate = dt * verhulst_rate(N[i-1],alpha,K)
        N[i] = N[i-1] + rate
    return N

#2nd order runga kutta
def runga_kutta_2(Nini,tini,alpha,K,dt,n):
    N = np.zeros((n))
    N[0] = Nini
    for i in xrange(1,n):
        k1 = dt * verhulst_rate(N[i-1],alpha,K)
        k2 = dt * verhulst_rate((N[i-1]+k1/2.0),alpha,K)
        N[i] = N[i-1] + k2
    return N

#4th order runga kutta
def runga_kutta_4(Nini,tini,alpha,K,dt,n):
    N = np.zeros((n))
    N[0] = Nini
    for i in xrange(1,n):
        k1 = dt * verhulst_rate(N[i-1],alpha,K)
        k2 = dt * verhulst_rate((N[i-1]+k1/2.0),alpha,K)
        k3 = dt * verhulst_rate((N[i-1]+k2/2.0),alpha,K)
        k4 = dt * verhulst_rate((N[i-1]+k3),alpha,K)
        N[i] = N[i-1] + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0
    return N

#euler
N_numerical[:,0] = euler(N_initial,0,alpha,K,dt,n)

#2nd order runga kutta    
N_numerical[:,1] = runga_kutta_2(N_initial,0,alpha,K,dt,n)

#4th order runga kutta    
N_numerical[:,2] = runga_kutta_4(N_initial,0,alpha,K,dt,n)
    
#analytical soln
N_analytical = verhulst_analytical(alpha,K,N_initial,t)

plt.figure(1,facecolor = 'white')
plt.plot(t,N_numerical[:,0],label='Euler Method')
plt.plot(t,N_numerical[:,1],label=r'$2^{nd}$ Order Runga Kutta')
plt.plot(t,N_numerical[:,2],label=r'$4^{th}$ Order Runga Kutta')
plt.plot(t,N_analytical,label='Analytical Solution',ls='--',lw = 2)
plt.legend(loc='best')
plt.xlim(0,250)
plt.xlabel('time [T]')
plt.ylabel('populations [# of individuals]')
plt.savefig('Euler_and_2_4Runga.png')
plt.show()

Problem 2

<br data-mce-bogus="1">
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import random

#Verhulst parameters
alpha = .1 #growth rate
K = 1000. #carrying capacity

#initial condition
N_initial = 1.0
runs = 100
num = int(K)-int(N_initial)+1

N = np.linspace(N_initial,K,num)
T = np.zeros((num,runs))
rnd = np.random.rand(num,runs)

T_avg = np.zeros((int(K)-int(N_initial)+1))

plt.figure(1,facecolor='white')

for i in xrange(0,runs):
    rate = alpha * N[1:] * (1.0 - N[1:] / K)
    dt = - np.log(1.0 - rnd[1:,i]) / rate
    T[1:,i] = np.cumsum(dt)
    T_avg += T[:,i]
    plt.plot(T,N,color ='b',alpha=0.1,lw = .5)
    
T_avg = T_avg / float(runs)
plt.plot(T_avg,N,color ='k',label ='Averaged Solution',lw =2)
N_analytical = (K * N_initial * np.exp(alpha*T_avg))/(K+N_initial*(np.exp(alpha*T_avg)-1.))    
plt.plot(T_avg,N_analytical,label='Analytical Solution',ls='--',lw = 2,color ='r')
plt.xlabel('t [days]')
plt.ylabel('N [population]')
plt.savefig('random.png')
plt.show()