Week 8: Solving Differentital Equations

First lets take a simple relation ship, $$ \frac{dx}{dt} = -\alpha x $$

To solve this all we would have to do is realize that this is exponential decay and the solution is $$x(t) = e^{-\alpha t} $$

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'
In [2]:
X=np.arange(0,10,0.01)

plt.plot(X,np.exp(-X))

plt.show()
In [5]:
def xfunc(x):
    return -x

We want to build a step in the derivative to approximate the exact solution

In [11]:
X=np.arange(0,2,0.01)

stepsize =0.4

x,y0=min(X),1

XX = [x]
YY = [y0]

while x< max(X):
    y_new = xfunc(YY[-1])*stepsize + YY[-1]
    YY.append(y_new)
    x+=stepsize
    XX.append(x)



plt.plot(X,np.exp(-X))
plt.plot(XX,YY,'.-','r')

plt.show()
In [15]:
X=np.arange(0,2,0.01)

for stepsize in np.linspace(0.05,0.5,5):
    
    x,y0=min(X),1
    XX = [x]
    YY = [y0]

    while x< max(X):
        y_new = xfunc(YY[-1])*stepsize + YY[-1]
        YY.append(y_new)
        x+=stepsize
        XX.append(x)
    plt.plot(XX,YY,'.-',label = str(round(stepsize,2)), alpha = 0.3)
    


plt.plot(X,np.exp(-X), label='exact solution')

plt.legend()

plt.show()

Lets Generalize this approach

$$\frac{dx}{dt}=f(x,t)$$

We want $f(x,t)$ to be dependent on position and time.

We want to plot a solution of this first order differential equation

We have the following inputs

  • we have a function $f(x,t)$
  • inital position at time zero
  • a stepsize
  • End time
    • In [16]:
    def xfunc(x,t):
    return -x

def ode_1d_solve(xfunc,X0,dt,T):
    Num = int(T//dt)
    x = np.ones(Num +1)*X0
    t = np.linspace(0,Num*dt,len(x))
    for count in range(Num):
        x[count+1]=x[count]+xfunc(x[count],t[count])*dt
        
    return t,x
    In [18]:
    t,x = ode_1d_solve(xfunc,1,0.2,2)

plt.plot(t,x,'.-')
plt.show()
    In [35]:
    %%time
X=np.arange(0,2,0.01)

t,x = ode_1d_solve(xfunc,1,0.02,2)

plt.plot(t,x,'.','r')


plt.plot(X,np.exp(-X), 'b',label='exact solution')

plt.legend()

plt.show()

    CPU times: user 200 ms, sys: 15.4 ms, total: 216 ms
Wall time: 242 ms
    In [33]:
    len(t)
    Out[33]:
    100000
    In [42]:
    %%time

def xfunc(x,t):
    return (np.exp(t/3)-x)/2



t,x = ode_1d_solve(xfunc,1,0.02,20)

plt.plot(t,x,'r','.')

plt.show()

    CPU times: user 169 ms, sys: 13.1 ms, total: 182 ms
Wall time: 198 ms

    Lets look at the logistic growth, in the continuous model

    $$ \frac{dx}{dt}=rx(K-x)$$
    In [46]:
    def logistic_func(x,t):
    return 0.15*x*(10-x)
    In [48]:
    t,x = ode_1d_solve(logistic_func,1,0.01,5)

plt.plot(t,x,'r')

plt.show()
    In [52]:
    plt.figure(figsize=(13,5))

for s in np.linspace(1,21,11):
    t,x = ode_1d_solve(logistic_func,s,0.001,6)
    
    plt.plot(t,x, label='starting value = '+str(round(s,1)))
    
plt.legend(loc='center left', bbox_to_anchor=(1,0.5))
plt.show()

    Lets look at 2-dimensional models

    Generalized preditor-prey model

    x being the prey, y being predator

    $$\frac{dx}{dt}=xf(x,y)$$

    $f(x,y)$ this is going to a factor that essentially describes how efficient the prey are

    $$\frac{dy}{dt}=yg(x,y)$$

    $g(x,y)$ this is going to a factor that describes the effect of the prey on the predators, as in carry capacity

    Lotka-Volterra Model

    $$\frac{dx}{dt}=x(b-py)$$$$\frac{dy}{dt}=y(rx-d)$$

    where $b,p,r,d >0$

    In [2]:
    b,p=1,1
r,d=1,2

def xfunc(x,y):
    return x*(b-p*y)
def yfunc(x,y):
    return y*(r*x-d)
    In [38]:
    T = 4
dt=0.001

X0,Y0 =1,0.5

Num = int(T//dt)

x = np.ones(Num +1)*X0
y = np.ones(Num +1)*Y0

t = np.linspace(0,Num*dt,len(x))

for count in range(Num):
    x[count+1]=x[count]+xfunc(x[count],y[count])*dt
    y[count+1]=y[count]+yfunc(x[count],y[count])*dt

    
plt.figure(figsize=(13,6))

plt.subplot(2,2,1)
plt.title('predator-prey dynamics')
plt.plot(x,y)
plt.plot(x[-1],y[-1],'x')
plt.plot(x[0],y[0],'o')

plt.ylabel('predator')

plt.subplot(2,2,2)
plt.plot(t,y)
plt.title('time-predator')
plt.plot(t[-1],y[-1],'x')
plt.plot(t[0],y[0],'o')
plt.ylabel('predator')

plt.subplot(2,2,3)
plt.plot(x,t)
plt.title('prey-time')
plt.plot(x[-1],t[-1],'x')
plt.plot(x[0],t[0],'o')
plt.xlabel('prey')

plt.show()
    In [42]:
    T = 30
dt=0.0001

X0,Y0 =1,0.5

Num = int(T//dt)

x = np.ones(Num +1)*X0
y = np.ones(Num +1)*Y0

t = np.linspace(0,Num*dt,len(x))

for count in range(Num):
    x[count+1]=x[count]+xfunc(x[count],y[count])*dt
    y[count+1]=y[count]+yfunc(x[count],y[count])*dt

    
plt.figure(figsize=(13,6))

plt.subplot(2,2,1)
plt.title('predator-prey dynamics')
plt.plot(x,y)
plt.plot(x[-1],y[-1],'x')
plt.plot(x[0],y[0],'o')

plt.ylabel('predator')

plt.subplot(2,2,2)
plt.plot(t,y)
plt.title('time-predator')
plt.plot(t[-1],y[-1],'x')
plt.plot(t[0],y[0],'o')
plt.ylabel('predator')

plt.subplot(2,2,3)
plt.plot(x,t)
plt.title('prey-time')
plt.plot(x[-1],t[-1],'x')
plt.plot(x[0],t[0],'o')
plt.xlabel('prey')

plt.show()
    In [5]:
    def Ode_2d_solve(xfunc,yfunc,x0,y0,dt,T):
    
    Num = int(T//dt)

    x = np.ones(Num +1)*x0
    y = np.ones(Num +1)*y0

    t = np.linspace(0,Num*dt,len(x))

    for count in range(Num):
        x[count+1]=x[count]+xfunc(x[count],y[count])*dt
        y[count+1]=y[count]+yfunc(x[count],y[count])*dt
        
    return t,x,y
    In [46]:
    b,p=1,1
r,d=1,2

t,x,y = Ode_2d_solve(xfunc,yfunc,1,0.3,0.00001,40)

plt.plot(x,y)
plt.title('predator-prey dynamics')
plt.xlabel('prey')
plt.ylabel('predator')
plt.show()
    In [51]:
    b,p=1,1
r,d=1,2

t,x,y = Ode_2d_solve(xfunc,yfunc,1,0.3,0.001,40)

plt.plot(x,y)
plt.title('predator-prey dynamics')
plt.xlabel('prey')
plt.ylabel('predator')
plt.show()
plt.figure(figsize=(13,6))
plt.plot(t,x,c='b',label='prey')
plt.plot(t,y,c='r',label='predator')
plt.legend()
plt.show()

    What if there was a limit to the food?

    Carrying Capacity on the prey!!

    We have to change the $\frac{dx}{dt}$

    $$\frac{dx}{dt}=x(b-py) =xb-pxy$$$$\frac{dx}{dt}= xb(c-x) - pxy$$$$\frac{dx}{dt}= bx-cx^2 - pxy$$
    In [3]:
    b,p,c=2,1,1
r,d = 5,2

def xfunc(x,y):
    return b*x-c*x*x-p*x*y
def yfunc(x,y):
    return y*(r*x-d)
    In [54]:
    t,x,y = Ode_2d_solve(xfunc,yfunc,1,0.3,0.0001,40)

plt.plot(x,y)
plt.title('predator-prey dynamics')
plt.xlabel('prey')
plt.ylabel('predator')
plt.show()
plt.figure(figsize=(13,6))
plt.plot(t,x,c='b',label='prey')
plt.plot(t,y,c='r',label='predator')
plt.legend()
plt.show()
    In [57]:
    x[5], y[5],xfunc(x[5],y[5])
    Out[57]:
    (1.0003498889322198, 0.30045037519040574, 0.6994443781263688)
    In [69]:
    T = 30
dt=0.0001

X0,Y0 =1,0.5

Num = int(T//dt)

noise_level = 8

x = np.ones(Num +1)*X0
y = np.ones(Num +1)*Y0

t = np.linspace(0,Num*dt,len(x))

for count in range(Num):
    x[count+1]=x[count]+np.random.normal(1,noise_level)*xfunc(x[count],y[count])*dt
    y[count+1]=y[count]+np.random.normal(1,noise_level)*yfunc(x[count],y[count])*dt

    
plt.figure(figsize=(13,6))

plt.subplot(2,2,1)
plt.title('predator-prey dynamics')
plt.plot(x,y)
plt.plot(x[-1],y[-1],'x')
plt.plot(x[0],y[0],'o')

plt.ylabel('predator')

plt.subplot(2,2,2)
plt.plot(t,y)
plt.title('time-predator')
plt.plot(t[-1],y[-1],'x')
plt.plot(t[0],y[0],'o')
plt.ylabel('predator')

plt.subplot(2,2,3)
plt.plot(x,t)
plt.title('prey-time')
plt.plot(x[-1],t[-1],'x')
plt.plot(x[0],t[0],'o')
plt.xlabel('prey')

plt.show()
    In [60]:
    YY=[]
for count in range(10000):
    YY.append(np.random.normal(1,noise_level))
    In [65]:
    plt.hist(YY,bins=100)
plt.show()
    In [70]:
    def Ode_2d_noise_solve(xfunc,yfunc,x0,y0,dt,T,noise_level):
    Num = int(T//dt)

    x = np.ones(Num +1)*x0
    y = np.ones(Num +1)*y0

    t = np.linspace(0,Num*dt,len(x))

    for count in range(Num):
        
        x[count+1]=x[count]+np.random.normal(1,noise_level)*xfunc(x[count],y[count])*dt
        y[count+1]=y[count]+np.random.normal(1,noise_level)*yfunc(x[count],y[count])*dt
        
    return t,x,y
    In [72]:
    t,x,y = Ode_2d_noise_solve(xfunc,yfunc,1,3,0.001,20,8)

plt.plot(x,y)
plt.title('predator-prey dynamics')
plt.xlabel('prey')
plt.ylabel('predator')
plt.show()
plt.figure(figsize=(13,6))
plt.plot(t,x,c='b',label='prey')
plt.plot(t,y,c='r',label='predator')
plt.legend()
plt.show()
    In [78]:
    t,x,y = Ode_2d_solve(xfunc,yfunc,1,3,0.0001,30)

t1,x1,y1 = Ode_2d_solve(xfunc,yfunc,1,2,0.0001,30)
t2,x2,y2 = Ode_2d_solve(xfunc,yfunc,1,1,0.0001,30)
plt.plot(x,y)
plt.plot(x1,y1)
plt.plot(x2,y2)
plt.title('predator-prey dynamics')
plt.xlabel('prey')
plt.ylabel('predator')
plt.show()


plt.figure(figsize=(13,6))
plt.plot(t,x,c='b',label='prey')
plt.plot(t,y,c='r',label='predator')
plt.plot(t,x1,c='b')
plt.plot(t,y1,c='r')
plt.plot(t,x2,c='b')
plt.plot(t,y2,c='r')

plt.legend()
plt.show()
    In [15]:
    b,p,c=2,1,1
r,d = 5,2

def xfunc(x,y):
    return b*x-c*x*x-p*x*y
def yfunc(x,y):
    return y*(r*x-d)



X = np.linspace(0.1,1,4)
Y = np.linspace(0.1,4,4)

X,Y = np.meshgrid(X,Y)
X = np.ravel(X)
Y = np.ravel(Y)
timelength = 20
stepsize =0.01

for count in range(len(X)):
    t1,x1,y1 = Ode_2d_solve(xfunc,yfunc,X[count],Y[count],stepsize,timelength)
    plt.plot(x1,y1,c='b',alpha=0.3)
    In [20]:
    X = np.linspace(0.1,1,4)
Y = np.linspace(0.1,4,4)

X,Y = np.meshgrid(X,Y)
X = np.ravel(X)
Y = np.ravel(Y)
timelength = 20
stepsize =0.01
plt.figure(figsize=(13,5))
for count in range(len(X)):
    t1,x1,y1 = Ode_2d_solve(xfunc,yfunc,X[count],Y[count],stepsize,timelength)
    plt.plot(t1,y1,c='b',alpha=0.3)
    plt.plot(t1,x1,c='r',alpha=0.2)
    In [23]:
    b,p=1,1
r,d=1,2

def xfunc(x,y):
    return x*(b-p*y)
def yfunc(x,y):
    return y*(r*x-d)



X = np.linspace(0.1,1,4)
Y = np.linspace(0.1,4,4)

X,Y = np.meshgrid(X,Y)
X = np.ravel(X)
Y = np.ravel(Y)
timelength = 20
stepsize =0.0001

for count in range(len(X)):
    t1,x1,y1 = Ode_2d_solve(xfunc,yfunc,X[count],Y[count],stepsize,timelength)
    plt.plot(x1,y1,c='b',alpha=0.3)
    In [24]:
    b,p,c=2,1,1
r,d = 5,2

def xfunc(x,y):
    return b*x-c*x*x-p*x*y
def yfunc(x,y):
    return y*(r*x-d)
    In [31]:
    for r  in range(10,0,-1):
    X = np.linspace(0.1,1,4)
    Y = np.linspace(0.1,4,4)

    X,Y = np.meshgrid(X,Y)
    X = np.ravel(X)
    Y = np.ravel(Y)
    timelength = 30
    stepsize =0.01

    for count in range(len(X)):
        t1,x1,y1 = Ode_2d_solve(xfunc,yfunc,X[count],Y[count],stepsize,timelength)
        plt.plot(x1,y1,c='b',alpha=0.3)
    plt.title('The r value is '+str(r))
    plt.show()
    In [ ]: