Introduction to Differential Equations
John Travis
Mississippi College
This worksheet investigates the nature of solutions to increasingly complex population growth models. Students can experiment with various parameters to find how each affects the nature of solutions.
t,p = var('t p')
a = 0 # Starting point for solution curve
b = 3 # Ending point for solution curve
p0 = 2 # Initial condition, p(a) = p0
f = -p0*p # Solving the first order DE P' = f
pmin = -3 # You may have to adjust the max and min values here to fit the window you want to see
pmax = 3
G = Graphics() # Set up a place to accumulate graphical objects
G = plot_slope_field(f, (t,a,b), (p,pmin,pmax)) # Note, this only involves the rhs f for P'=f
P = function('P')(t) # Set up a symbolic object P for the exact solution found below
DE = diff(P, t, 1) - (-p0*P)
soln = desolve( DE, P, ics=[a,p0] )
G += plot(soln,(t,a,b))
show(G)
show(DE," = 0")
show("P(t) = ",soln)
|
# Now, let's plot a number of solution curves for various initial starting points
t,p = var('t p')
a = 0 # Starting point for solution curve
b = 3 # Ending point for solution curve
k = -2 # Growth parameter if k>0, Decay parameter if k<0
f = k*p # Solving the first order DE P' = f
pmin = -3 # You may have to adjust the max and min values here to fit the window you want to see
pmax = 3
G = Graphics() # Set up a place to accumulate graphical objects
G = plot_slope_field(f, (t,a,b), (p,pmin,pmax)) # Note, this only involves the rhs f for P'=f
P = function('P')(t) # Set up a symbolic object P for the exact solution found below
for j in range(10): # Do 10 solution curves for 10 equally spaced starting values
delp = (pmax-pmin)/9
p0 = pmin+j*delp
DE = diff(P, t, 1) - (k*P)
soln = desolve( DE, P, ics=[a,p0] )
G += plot(soln,(t,a,b))
show(G)
show(DE," = 0")
show("P(t) = ",soln)
|
# Now, let's make this graphic interactive based upon the parameter k
t,p = var('t p')
a = 0 # Starting point for solution curve
b = 3 # Ending point for solution curve
@interact
def _(k=slider(-3,1,0.1,-2)):
f = k*p # Solving the first order DE P' = f
pmin = -3 # You may have to adjust the max and min values here to fit the window you want to see
pmax = 3
G = Graphics() # Set up a place to accumulate graphical objects
G = plot_slope_field(f, (t,a,b), (p,pmin-1,pmax+1)) # Note, this only involves the rhs f for P'=f
P = function('P')(t) # Set up a symbolic object P for the exact solution found below
for j in range(10): # Do 10 solution curves for 10 equally spaced starting values
delp = (pmax-pmin)/9
p0 = pmin+j*delp
DE = diff(P, t, 1) - (k*P)
soln = desolve( DE, P, ics=[a,p0] )
G += plot(soln,(t,a,b))
show(G)
show(DE," = 0")
show("P(t) = ",soln)
Click to the left again to hide and once more to show the dynamic interactive window |
# Let's change the Differential Equation to make it "First order linear"
t,p = var('t p')
a = 0 # Starting point for solution curve
b = 3 # Ending point for solution curve
@interact
def _(p0 = input_box(default=1,width=10)): # Initial condition, p(a) = p0
f = t-p0*p # Solving the first order DE P' = f
pmin = -3 # You may have to adjust the max and min values here to fit the window you want to see
pmax = 3
G = Graphics() # Set up a place to accumulate graphical objects
G = plot_slope_field(f, (t,a,b), (p,pmin,pmax)) # Note, this only involves the rhs f for P'=f
P = function('P')(t) # Set up a symbolic object P for the exact solution found below
DE = diff(P, t, 1) - (t-p0*P)
soln = desolve( DE, P, ics=[a,p0] )
G += plot(soln,(t,a,b),color='red', thickness=5, alpha=0.5) + points((a,p0),color='green',size=50)
show(G)
show(DE," = 0")
show("P(t) = ",soln)
Click to the left again to hide and once more to show the dynamic interactive window |
# John Travis
# Mississippi College
# here is a numerical solver using splines on numerical trajectories.
t,P = var('t,P') # define an independent variable t
var('k,M,N,C')
models = ['Exponential Growth','Logistic Growth','Logistic Growth with Sustainability','Logistic Growth with Harvesting']
@interact
def _(k=slider(0.1,1.0,0.1,0.5,label='$k$'),
model=models,
M=input_box(1,width=10,label='Carrying capacity'),
N=input_box(4/10,width=10,label='Sustainability min.'),
C=input_box(1/10,width=10,label='Harvesting level')):
t1 = 5
if (model==models[0]):
method=0
elif (model==models[1]):
method=1
elif (model==models[2]):
method=2
else:
method=3
F = [k*P,k*P*(1-P/M),k*P*(1-P/M)*(P-N),k*P*(1-P/M)-C] # Here is the complicated DE P'=f
flatex = ['$dP/dt = kP$','$dP/dt = kP(1-P/M)$','$dP/dt = kP(1-P/M)(P-N)$','$dP/dt = kP(1-P/M)-C$']
v = []
pts = Graphics()
# G = plot_slope_field(f, (t,a,b), (p,pmin,pmax)) # Note, this only involves the rhs f for P'=f
P0 = [n/10 for n in range(1, 21)] # initial population choices
f = F[method]
for p0 in P0:
soln = desolve_rk4(f,P,ics=[0,p0],ivar=t,end_points=t1,step=0.25)
pts += plot(spline(soln),0,t1)
v.append(pts)
pts.show(ymin=0,ymax=4)
# animate(v,xmin=0,xmax=t1,ymin=0,ymax=2*M).show()
pretty_print(html('<center>Solution trajectories for\n'+flatex[method]+'\nusing various initial conditions.</center>'))
Click to the left again to hide and once more to show the dynamic interactive window |
||||||||||||||||||||||||||
# John Travis
# Mississippi College
# here is a numerical solver using splines on numerical trajectories.
t,P = var('t,P') # define an independent variable t
var('k,M,N,C')
models = ['Exponential Growth','Logistic Growth','Logistic Growth with Sustainability','Logistic Growth with Harvesting']
@interact
def _(k=slider(0.1,1.0,0.1,0.5,label='$k$'),
model=models,
M=input_box(1,width=10,label='Carrying capacity'),
N=input_box(4/10,width=10,label='Sustainability min.'),
C=input_box(1/10,width=10,label='Harvesting level')):
t1 = 5
if (model==models[0]):
method=0
elif (model==models[1]):
method=1
elif (model==models[2]):
method=2
else:
method=3
F = [k*P,k*P*(1-P/M),k*P*(1-P/M)*(P-N),k*P*(1-P/M)-C] # Here is the complicated DE P'=f
flatex = ['$dP/dt = kP$','$dP/dt = kP(1-P/M)$','$dP/dt = kP(1-P/M)(P-N)$','$dP/dt = kP(1-P/M)-C$']
v = []
pts = Graphics()
P0 = [n/10 for n in range(1, 21)] # initial population choices
f = F[method]
for p0 in P0:
soln = desolve_rk4(f,P,ics=[0,p0],ivar=t,end_points=t1,step=0.25)
pts += plot(spline(soln),0,t1)
v.append(pts)
pts.show(ymin=0,ymax=4)
# animate(v,xmin=0,xmax=t1,ymin=0,ymax=2*M).show()
pretty_print(html('<center>Solution trajectories for\n'+flatex[method]+'\nusing various initial conditions.</center>'))
Click to the left again to hide and once more to show the dynamic interactive window |