Math 1003, Week 8, Lectures 15 and 16

Lecturer: Eduardo G. Altmann, http://www.maths.usyd.edu.au/u/ega/

1. Application of Separable Equations

1.2 Radioactive decay

$$ \frac{d N}{d t} = k N $$

import matplotlib.pyplot as plt
from scipy import *
from scipy import integrate
from scipy.integrate import ode
import numpy as np

%load_ext autoreload
%autoreload 2

%matplotlib inline

fig = plt.figure(num=1,figsize=(10,10))
ax=fig.add_subplot(111)

xmin=0;xmax=100;ymin=0;ymax=100;
gridx=30;gridy=30;
x=np.array(range(100))/100.;x=xmin+(xmax-xmin)*x;


plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel(r"$t$",fontsize=18)
plt.ylabel(r"$N$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

k=np.log(2)/30.0

#Vector field
X,Y = np.meshgrid( np.linspace(xmin,xmax,gridx),np.linspace(ymin,ymax,gridy) )
dX = 1
dY = -k*Y

#Normalize arrows
N = np.sqrt(dX**2+dY**2)  
U, V = dX/N, dY/N
ax.quiver( X,Y,U, V,minshaft=5,width=0.005,headwidth=0,color="black")


# Ploting solutions derived in classroom
f=100*exp(-k*x)
ax.plot(x,f)

f=50*exp(-k*x)
ax.plot(x,f)

plt.show()

1.2 Physiological Data

$$ \frac{d R}{dS} = k R/S$$

fig = plt.figure(num=1,figsize=(10,10))
ax=fig.add_subplot(111)

xmin=0.01;xmax=10;ymin=0;ymax=10;
gridx=30;gridy=30;
x=np.array(range(100))/100.;x=xmin+(xmax-xmin)*x;


plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel(r"$S$",fontsize=18)
plt.ylabel(r"$R$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

k=2

#Vector field
X,Y = np.meshgrid( np.linspace(xmin,xmax,gridx),np.linspace(ymin,ymax,gridy) )
dX = 1
dY = k*Y/X

#Normalize arrows
N = np.sqrt(dX**2+dY**2)  
U, V = dX/N, dY/N
ax.quiver( X,Y,U, V,minshaft=5,width=0.005,headwidth=0,color="black")


# Ploting solutions derived in classroom
f=2*x**k
ax.plot(x,f)

f=0.1*x**k
ax.plot(x,f)

plt.show()

1.3 Classical Mechanics

$$\frac{dv}{dt} = g - kv $$

fig = plt.figure(num=1,figsize=(10,10))
ax=fig.add_subplot(111)

xmin=0.0;xmax=4;ymin=0;ymax=10;
gridx=30;gridy=30;
x=np.array(range(100))/100.;x=xmin+(xmax-xmin)*x;


plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel(r"$S$",fontsize=18)
plt.ylabel(r"$R$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

k=2
g=9.8

#Vector field
X,Y = np.meshgrid( np.linspace(xmin,xmax,gridx),np.linspace(ymin,ymax,gridy) )
dX = 1
dY = 9.8-2*Y

#Normalize arrows
N = np.sqrt(dX**2+dY**2)  
U, V = dX/N, dY/N
ax.quiver( X,Y,U, V,minshaft=5,width=0.005,headwidth=0,color="black")

# Ploting solutions derived in classroom
v0=0
f=(g/k)*(1-(1-k*v0/g)*(np.exp(-k*x)))
ax.plot(x,f)

v0=8
f=(g/k)*(1-(1-k*v0/g)*(np.exp(-k*x)))
ax.plot(x,f)

plt.show()

2 Models of Growth

2.1 Linear vs. Exponential growth

$$\frac{dy}{dt} = a \Rightarrow y(t)= at+y(0) $$

vs.

$$\frac{dy}{dt} = ky \Rightarrow y(0) e^{k t}$$

fig=plt.figure(figsize=(20,10))

xmin=0.0;xmax=10;ymin=0.1;ymax=10;
x=np.array(range(1000))/1000.;x=xmin+(xmax-xmin)*x;

plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)
#plt.yscale("log")
plt.xlabel(r"$t$",fontsize=30)
plt.ylabel(r"$y$",fontsize=30)

plt.tick_params(axis='both', which='major', labelsize=30)

plt.plot(x,1*x+0,linewidth=2,label="Linear a=1, y(0)=0")
plt.plot(x,0.5*x+2,linewidth=2,label="Linear a=1/2, y(0)=2")
plt.plot(x,exp(1*x),'--',linewidth=2,label="Exponential k=1, y(0)=1")
plt.plot(x,0.1*exp(0.1*x),'--',linewidth=2,label="Exponential k=0.1, y(0)=0.1")
plt.legend(loc=0,fontsize=30)
plt.show(fig)

########################################################################################
fig=plt.figure(figsize=(20,10))
#xscvale("log")
xmin=0.0;xmax=100;ymin=0.1;ymax=100;
x=np.array(range(1000))/1000.;x=xmin+(xmax-xmin)*x;

plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)
plt.xlabel(r"$t$",fontsize=30)
plt.ylabel(r"$y$",fontsize=30)

plt.tick_params(axis='both', which='major', labelsize=30)

plt.plot(x,1*x+0,linewidth=2,label="Linear a=1, y(0)=0")
plt.plot(x,0.5*x+2,linewidth=2,label="Linear a=1/2, y(0)=2")
plt.plot(x,exp(1*x),'--',linewidth=2,label="Exponential k=1, y(0)=1")
plt.plot(x,0.1*exp(0.1*x),'--',linewidth=2,label="Exponential k=0.1, y(0)=0.1")
#plt.legend(loc=0,fontsize=30)
plt.show(fig)

#################################################################################

fig=plt.figure(figsize=(20,10))
#xscvale("log")
xmin=0.0;xmax=100;ymin=0.1;ymax=10000;
x=np.array(range(1000))/1000.;x=xmin+(xmax-xmin)*x;

plt.xlim(xmin,xmax)
plt.ylim(ymin,ymax)
plt.yscale("log")
plt.xlabel(r"$t$",fontsize=30)
plt.ylabel(r"$y$",fontsize=30)

plt.tick_params(axis='both', which='major', labelsize=30)

plt.plot(x,1*x+0,linewidth=2,label="Linear a=1, y(0)=0")
plt.plot(x,0.5*x+2,linewidth=2,label="Linear a=1/2, y(0)=2")
plt.plot(x,exp(1*x),'--',linewidth=2,label="Exponential k=1, y(0)=1")
plt.plot(x,0.1*exp(0.1*x),'--',linewidth=2,label="Exponential k=0.1, y(0)=0.1")
#plt.legend(loc=0,fontsize=30)
plt.show(fig)

2.2 Logistic growth

$$ \frac{dx}{dt} = a x (b-x) $$

fig = plt.figure(num=1,figsize=(10,10))
ax=fig.add_subplot(111)

xmin=0.0;xmax=10;ymin=0;ymax=2;
gridx=30;gridy=30;
x=np.array(range(100))/100.;x=xmin+(xmax-xmin)*x;


plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel(r"$t$",fontsize=18)
plt.ylabel(r"$x$",fontsize=18)
plt.tick_params(axis='both', which='major', labelsize=16)

#ab=k
#K=x(0/(b-x(0)))
a=1
b=1
k=1

#Vector field
X,Y = np.meshgrid( np.linspace(xmin,xmax,gridx),np.linspace(ymin,ymax,gridy) )
dX = 1
dY = a*Y*(b-Y)

#Normalize arrows
N = np.sqrt(dX**2+dY**2)  
U, V = dX/N, dY/N
ax.quiver( X,Y,U, V,minshaft=5,width=0.005,headwidth=0,color="black")

# Ploting solutions derived in classroom
x0=1.5
C=(b-x0)/x0
f=b/(1+C*exp(-k*x))
ax.plot(x,f,linewidth=2)


x0=0.1
C=(b-x0)/x0
f=b/(1+C*exp(-k*x))
ax.plot(x,f,linewidth=2)

plt.show()