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>]
In [ ]: