352 - Topic 01 - Population Models

109 days ago by Professor352

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>')) 
       
$k$ 
model 
Carrying capacity 
Sustainability min. 
Harvesting level 

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