Week 9: Solving Non-linear Differentital Equations and 3 dimensions

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import Pythonic_files as pf
%config InlineBackend.figure_format='retina'

Nonlinear Pendulum

Using Newtons Second Law $F=ma$, we get

$$ -bl\frac{d\theta}{dt} -mg\sin\theta = ml \frac{d^2\theta}{dt^2}$$

where $b>0$, where $b$ is damping constant

Usually this, but we are getting funky $$ \frac{g}{l}\sin\theta = \frac{d^2\theta}{dt^2}$$

Change the variable to make things easy, lets just say that $v$ is the speed of the pendulum in angular

$v= \frac{d\theta}{dt}$, If we differentiate this function $ \frac{dv}{dt} = \frac{d^2\theta}{dt^2} $

We now have two differential equations, $$ \frac{d\theta}{dt}=v$$

$$ \frac{dv}{dt} = \frac{d^2\theta}{dt^2} = -\frac{b}{m}v-\frac{g}{l}sin\theta$$
In [2]:
b=2
g=9.81
l=1
m=3

def xfunc(x,y):
    return y

def yfunc(x,y):
    return -b/m*y-g/l*np.sin(x)
In [7]:
t,x,y = pf.Ode_2D_solve(xfunc,yfunc,2,1,0.01,10)

plt.plot(x,y,c='b')
plt.plot(x[0],y[0],'x',c='g')
plt.plot(x[-1],y[-1],'o',c='r')
Out[7]:
[<matplotlib.lines.Line2D at 0x7feb98577d30>]
In [33]:
plt.figure(figsize=(13,8))

b=3
g=9.81
l=2
m=4

def xfunc(x,y):
    return y

def yfunc(x,y):
    return -b/m*y-g/l*np.sin(x)

X=np.linspace(-13,13,30)
Y=np.linspace(-2,2,9)

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)):
    t,x,y = pf.Ode_2D_solve(xfunc,yfunc,X[count],Y[count],stepsize,timelength)
    plt.plot(x,y,c='b',alpha=0.1)
    
plt.show()

Three dimensional models

Wolves, Moose and Tree rings on Isle Royal, Science 1994

  • $x(t)=$ Population of Trees at time t
  • $y(t)=$ Population of Moose at time t
  • $z(t)=$ Population of Wolves at time t
$$ \frac{dx}{dt} = x(1-x) -xy $$$$ \frac{dy}{dt} = y(1-y) +xy -yz $$$$ \frac{dz}{dt} = z(1-z) +yz $$
In [34]:
def xfunc(x,y,z):
    return x*(1-x)-x*y

def yfunc(x,y,z):
    return y*(1-y) +x*y -y*z

def zfunc(x,y,z):
    return z*(1-z) +y*z
In [37]:
def Ode_3D_solve(xfunc,yfunc,zfunc,x0,y0,z0,dt,T):
    
    Num = int(T//dt)


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

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

    
    for count in range(Num):
        x[count +1] = x[count]+xfunc(x[count],y[count],z[count])*dt
        y[count +1] = y[count]+yfunc(x[count],y[count],z[count])*dt
        z[count +1] = z[count]+zfunc(x[count],y[count],z[count])*dt
    return t,x,y,z
In [55]:
fig = plt.figure()
ax= plt.axes(projection = "3d")

t,x,y,z = Ode_3D_solve(xfunc,yfunc,zfunc,10,2,1,0.01,1)

ax.plot3D(x,y,z)
ax.scatter3D(x[0],y[0],z[0],c='g')
ax.scatter3D(x[-1],y[-1],z[-1],c='r')
ax.set_xlabel('Trees')
ax.set_ylabel('Moose')
ax.set_zlabel('Wolves')
ax.view_init(30,55)

print(x[-1],y[-1],z[-1])

0.6554022196797712 0.7413505915288129 2.1456449881026396
In [58]:
plt.figure(figsize=(13,6))
plt.plot(t,x,c='b',label='Trees')
plt.plot(t,y,c='r',label='Moose')
plt.plot(t,z,c='g',label='Wolves')
plt.legend()
plt.ylim(0,5)
plt.show()
In [60]:
def xfunc(x,y,z):
    return 2*x*(1-x)-x*y

def yfunc(x,y,z):
    return y*(1-y) +x*y -y*z

def zfunc(x,y,z):
    return z*(1-z) +y*z

fig = plt.figure()
ax= plt.axes(projection = "3d")

t,x,y,z = Ode_3D_solve(xfunc,yfunc,zfunc,10,2,1,0.01,5)

ax.plot3D(x,y,z)
ax.scatter3D(x[0],y[0],z[0],c='g')
ax.scatter3D(x[-1],y[-1],z[-1],c='r')
ax.set_xlabel('Trees')
ax.set_ylabel('Moose')
ax.set_zlabel('Wolves')
ax.view_init(30,55)

print(x[-1],y[-1],z[-1])

plt.figure(figsize=(13,6))
plt.plot(t,x,c='b',label='Trees')
plt.plot(t,y,c='r',label='Moose')
plt.plot(t,z,c='g',label='Wolves')
plt.legend()
plt.ylim(0,5)
plt.show()

0.8056465524746529 0.39442100973674393 1.3892551696391526

Lorentz Equation

$$ \frac{dx}{dt} = \sigma (y-x) $$$$ \frac{dy}{dt} = \rho x - y -xz $$$$ \frac{dz}{dt} = -\beta z +xy $$
In [76]:
s=10
r=28
b=8/3

def xfunc(x,y,z):
    return s*(y-x)

def yfunc(x,y,z):
    return r*x-y-x*z

def zfunc(x,y,z):
    return -b*z+x*y


fig = plt.figure(figsize=(13,13))
ax= plt.axes(projection = "3d")

t,x,y,z = Ode_3D_solve(xfunc,yfunc,zfunc,1,1,1,0.01,100)

ax.plot3D(x,y,z,linewidth=0.5)
ax.view_init(30,25)
In [78]:
for g in range(0,90,5):
    fig = plt.figure(figsize=(13,13))
    ax= plt.axes(projection = "3d")
    
    ax.plot3D(x,y,z,linewidth=0.5)
    ax.view_init(30,g)
    
    plt.show()
In [82]:
plt.figure(figsize=(13,5))
plt.plot(t,x)
Out[82]:
[<matplotlib.lines.Line2D at 0x7feb88f64190>]
In [83]:
t,x,y,z = Ode_3D_solve(xfunc,yfunc,zfunc,1,1,1,0.01,100)
t1,x1,y1,z1 = Ode_3D_solve(xfunc,yfunc,zfunc,1.1,1,1,0.01,100)

plt.figure(figsize=(13,5))
plt.plot(t,x,c='b')
plt.plot(t1,x1,c='r')
Out[83]:
[<matplotlib.lines.Line2D at 0x7feb60327b80>]

Disease modelling: the SIR model

We can split a population in to different "compartments"

  • S: the susceptible, the people can catch the disease
  • I: Infectious people, people who presently have the disease
  • R: People who are not susceptible or Infectious
  • N: the total population

Lets look at the differential equations

$$\frac{ds}{dt} = \alpha -\beta s i $$$$\frac{di}{dt} = \beta s i -\gamma i$$$$\frac{dr}{dt} = \gamma i -\alpha $$

The total population $=N$, when we normalize we had $N=S+I+R$ and then we get after normalization $1=s+i+r$

$$r=1-i-s$$$$\frac{dr}{dt} = -\frac{di}{dt} - \frac{ds}{dt} $$

Hence our model,

$$\frac{ds}{dt} = -\beta s i $$$$\frac{di}{dt} = \beta s i -\gamma i$$
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import Pythonic_files as pf
import seaborn as sns

%config InlineBackend.figure_format='retina'

For our first example lets use $\beta = 2$ and $ \gamma = \frac{1}{4}$

In [2]:
beta , gamma = 2,1/4

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i
In [3]:
I0=0.001
S0=1-I0

t,s,i = pf.Ode_2D_solve(sfunc,ifunc,S0,I0,0.01,30)
In [4]:
plt.plot(t,s,label='Susceptible')
plt.plot(t,i,label='Infectious')
plt.plot(t,1-s-i,label='Recovered')


plt.legend()
plt.show()
In [37]:
plt.figure(figsize=(13,6))

beta , gamma = 0.7,1/14

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i

I0=0.001
S0=1-I0
dt=0.01

t,s,i = pf.Ode_2D_solve(sfunc,ifunc,S0,I0,dt,60)

plt.plot(t,s,label='Susceptible')
plt.plot(t,i,label='Infectious')
plt.plot(t,1-s-i,label='Recovered')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

plt.legend()
plt.show()
In [27]:
plt.figure(figsize=(13,6))
plt.subplot(2,2,1)
plt.plot(s,i)
plt.xlabel('susceptible')
plt.ylabel('infectious')
plt.plot(s[0],i[0],'o',c='g')
plt.plot(s[-1],i[-1],'x',c='r')

plt.subplot(2,2,2)
plt.plot(t,i)
plt.xlabel('time')
plt.ylabel('infectious')
plt.plot(t[0],i[0],'o',c='g')
plt.plot(t[-1],i[-1],'x',c='r')

plt.subplot(2,2,3)
plt.plot(s,t)
plt.xlabel('time')
plt.ylabel('susceptible')
plt.plot(s[0],t[0],'o',c='g')
plt.plot(s[-1],t[-1],'x',c='r')

plt.show()

How do you figure out the $\beta$ from data

We know $\frac{ds}{dt} = -\beta s i $, thats the infinitesimal derivative, so an approximation is $\frac{\Delta s}{\Delta t} = -\beta si$

Or

$$ \beta= - \frac{\frac{\Delta s}{\Delta t}}{si}$$

We can use for each data point a $\beta$ representation, so a collection of them will be a distribution of the $\beta$'s

In [53]:
Days=[]
for x in range(len(t)):
    if t[x]//1-t[x-1]//1>=1:
        Days.append(x)
In [57]:
Betas = []
for x in Days[1:]:
    b = -(s[x]-s[x-1])/dt/(s[x]*i[x])
    Betas.append(b)
In [61]:
Betas

plt.hist(Betas)
plt.show()
In [71]:
plt.figure(figsize=(13,6))

beta , gamma = 0.3365,1/20

def sfunc(s,i):
    return -beta*s*i
def ifunc(s,i):
    return beta*s*i-gamma*i

I0=0.001
S0=1-I0
dt=0.01

t,s,i = pf.Ode_2D_solve(sfunc,ifunc,S0,I0,dt,100)

plt.plot(t,s,label='Susceptible')
plt.plot(t,i,label='Infectious')
plt.plot(t,1-s-i,label='Recovered')

plt.xlabel('time in days/months/mucci')
plt.ylabel('Proportion of Population')

plt.legend()
plt.show()
Days=[]
for x in range(len(t)):
    if t[x]//1-t[x-1]//1>=1:
        Days.append(x)
Betas = []
for x in Days[1:]:
    b = -(s[x]-s[x-1])/dt/(s[x]*i[x])
    Betas.append(b)
    
sns.distplot(Betas)
plt.show()
In [ ]: