#DMB14.txt: Maple Code for Dr. Z.'s Dynamical Models in Biology class, Lecture 14, Oct. 17, 2025 with(linalg): with(LinearAlgebra): print(`For a list of the Main procedures type: Help14(); for help with a specific procedure type: Help14(ProcedureName); for example Help14(Feig);`): Help14a:=proc() if nargs=0 then print(`The other procedures are: `): else Help14(args): fi: end: Help14:=proc() if nargs=0 then print(`The available procedures are: ChaosStory, Feig, HW, NicholsonBailey, NicholsonBaileyM, Orb, ORB, Orbk, RecToT, RecToTs, SS, SSg, SSr, SSS, SSSg, RR, RT `): elif nargs=1 and args[1]=ChaosStory then print(`ChaosStory(): This illustarates the period-doubling phenomenon. Try:`): print(`ChaosStory(); `): elif nargs=1 and args[1]=Feig then print(`Feig(): The list (from wikipedia) of length 8, of the period-doubling cut-offs`): print(`of the logistic equation x(n+1)=r*x(n)*(1-x(n)). Calling this list L, i.e.`): print(`L:=Feig(): The r below L[1] every orbit tends to ONE ultimate value`): print(`for r between L[1] and L[2], in the long run it alternates between 2 ultimate values`): print(`between L[2] and L[3] , in the long run it alternates between 4=2^2 ultimate values`): print(`Between L[i] and L[i+1], in the long run it alternates between 2^i ultimate values. Try:`): print(`Feig():`): elif nargs=1 and args[1]=HW then print(`HW(u,v): The Hardy-Weinberg unerlying transformation witu (u,v,w), Eqs. (53a,53b,53c) in Edelestein-Keshet Ch. 3 using the fact that u+v+w=1. Try:`): print(`HW(u,v);`): elif nargs=1 and args[1]=NicholsonBailey then print(`NicholsonBailey(L,a,c,N,P): The discrete-time, double-species dynamical model of Nicholson and Bailey (1935), given by Eqs. (21a)(21b) in Edelstein-Keshet section 3.2 (p. 81)`): print(`where the variables are N (hosts) and parasites (P) and the parameters are L (called Lambda there), a, and c`): print(`Try:`): print(`NicholsonBailey(L,a,c,N,P);`): print(`NicholsonBailey(2,0.068,1,N,P);`): elif nargs=1 and args[1]=NicholsonBaileyM then print(`NicholsonBaileyM(a,r,K,N,B): The discrete-time, double-species dynamical model of the MODIFIED Nicholson and Bailey model (1935), given by Eqs. (28a)(28b) in Edelstein-Keshet section 3.4 (p. 84)`): print(`where the variables are N (hosts) and parasites (P) and the parameters are r and K`): print(`Try:`): print(`NicholsonBaileyM(r,a,K,N,P);`): print(`NicholsonBaileyM(0.5,0.1,14,N,P);`): elif nargs=1 and args[1]=Orb then print(`Orb(f,x,x0,K1,K2): Given a function f of the variable x, an initial value x(0) finds the portion of the orbit`): print(`of the recurrence a(n+1)= f(a(n)), in other words [a(K1),a(K1+1), ..., a(K2)]. If you want the full orbit take K1=1.`): print(`For example, try:`): print(`Orb(5/2*x*(1-x),x,0.2,1000,1010);`): elif nargs=1 and args[1]=ORB then print(`ORB(F,x,x0,K1,K2): Inputs a transformation F in the list of variables x with initial point pt, outputs the trajectory`): print(`of the discrete dynamical system (i.e. solutions of the difference equation): x(n)=F(x(n-1)) with x(0)=x0 from n=K1 to n=K2.`): print(`For the full trajectory (from n=0 to n=K2), use K1=0. Try:`): print(`ORB([5/2*x*(1-x)],[x], [0.5], 1000,1010);`): print(`ORB([(1+x+y)/(2+x+y),(1+x)/(1+y)],[x,y], [2.,3.], 1000,1010);`): elif nargs=1 and args[1]=Orbk then print(`Orbk(k,z,f,INI,K1,K2): Given a positive integer k, a letter (symbol), z, an expression f of z[1], ..., z[k] (representing a multi-variable function of the variables z[1],...,z[k]`): print(`a vector INI representing the initial values [x[1],..., x[k]], and (in applications) positive integres K1 and K2, outputs the`): print(`values of the sequence starting at n=K1 and ending at n=K2. of the sequence satisfying the difference equation`): print(`x[n+k]=f(x[n+k-1],x[n+k-2],..., x[n]):`): print(`This is a generalization to higher-order difference equation of procedure Orb(f,x,x0,K1,K2). For example`): print(`Orbk(1,z,5/2*z[1]*(1-z[1]),[0.5],1000,1010); should be the same as`): print(`Orb(5/2*z[1]*(1-z[1]),z[1],[0,5],1000,1010);`): print(`Try: `): print(`Orbk(2,z,(5/4)*z[1]-(3/8)*z[2],[1,2],1000,1010);`): elif nargs=1 and args[1]=OrbkF then print(`OrbkF(k,z,f,INI,K1,K2): Like Orbk(k,z,f,INI,K1,K2) but in floating-point`): print(`OrbkF(1,z,5/2*z[1]*(1-z[1]),[0.5],1000,1010); should be the same as`): print(`OrbkF(5/2*z[1]*(1-z[1]),z[1],[0.5],1000,1010);`): print(`Try: `): print(`OrbkF(2,z,(5/4)*z[1]-(3/8)*z[2],[1,2],1000,1010);`): elif nargs=1 and args[1]=RecToT then print(`RecToT(k,z,f,INI): Given a k-th order recurrence a(n+k)=f(a(n+k-1), ..., a(n)) where f is written in terms of z[1], ..., z[k], and initial conditions`): print(`INI=[a(1),..., a(k)], outputs the transformation from F from R^k to R^k such that`): print(`if x(n)=[a(n+k-1), ..., a(n)] then x(n+1)= F(x(n)). Try:`): print(`RecToT(2,z,(z[1]+z[2])/(z[1]+3*z[2]),[1,2]);`): elif nargs=1 and args[1]=RecToTs then print(`RecToTs(k,z,f,INI): Given a k-th order recurrence a(n+k)=f(a(n+k-1), ..., a(n)) where f is written in terms of z[1], ..., z[k]`): print(` outputs the transformation from F from R^k to R^k such that`): print(`if x(n)=[a(n+k-1), ..., a(n)] then x(n+1)= F(x(n)). Try:`): print(`RecToTs(2,z,(z[1]+z[2])/(z[1]+3*z[2]));`): elif nargs=1 and args[1]=RR then print(`RR(var,K): A random rational function with numerator and denominator of degree 1 in the list of variables var, where the ocefficinets are random from 1 to K Try:`): print(`RR([x,y,z],10);`): print(`is for generating examples`): print(`Try:`): print(`RR([x,y],10);`): elif nargs=1 and args[1]=RT then print(`RT(var,K): A random rational transformation of numerator and denominator degrees 1 from R^k to R^k (where k=nops(var), with postiive integer coefficients from 1 to K The inpus are a list of variables x and a posi. integer K`): print(`is for generating examples`): print(`Try:`): print(`RT([x,y],10);`): elif nargs=1 and args[1]=SS then print(`SS(f,x): Given a function f of x outputs all the fixed points. Try:`): print(`SS(5/2*x*(1-x),x);`): elif nargs=1 and args[1]=SSg then print(`SSg(T,x): Given a transormation T in the list of variables, finds all the steady-states (in floating point). Try:`): print(`SSg([(1+x)/(3+y),(3+x)/(1+y)],[x,y]);`): elif nargs=1 and args[1]=SSr then print(`SSr(f,x): Given the underlying function of a recurrence a(n+k)=f(a(n+k-1), ..., a(n)) where f is expressed in terms of x[1],..., x[k] (and k:=nops(k)`): print(`Finds all the steady-states. Try:`): print(`SSr((3+x)/(1+y),[x,y]);`): elif nargs=1 and args[1]=SSS then print(`SSS(f,x): inputs a function f of x, outputs all the STABLE steady-states of the associated recurrence a(n+1)=f(a(n)). Try:`): print(`SSS(5/2*x*(1-x),x);`): elif nargs=1 and args[1]=SSSg then print(`SSSg(F,x): Given a transformation F in the list of variables finds all the STABLE fixed point of the transformation x->F(x), i.e. the set of solutions of`): print(`the system {x[1]=F[1], ..., x[k]=F[k]}} that are stable. Try:`): print(`SSSg([5/2*x*(1-x)],[x]]);`): print(`SSSg([(1+x+y)/(2+3*x+y), (3+x+2*y)/(5+x+3*y)],[x,y]);`): elif nargs=1 and args[1]=ToSysthen then print(`ToSys(k,z,f): converts the kth order difference equation a(n+k)=f(a[n+k-1],a[n+k-2],...a[n]) to a first-order system`): print(`where f is a function of z[1], z[2], ..., z[k]. Try:`): print(`ToSys(2,z,z[1]+z[2]);`): else print(`There is no Help for`, args): fi: end: #Orb(f,x,x0,K1,K2): Given a function f of the variable x, an initial value x(0) finds the portion of the orbit #of the recurrence a(n+1)= f(a(n)), in other words [a(K1),a(K1+1), ..., a(K2)]. If you want the full orbit take K1=1. #For example, try: #Orb(5/2*x*(1-x),x,0.2,1000,1010); Orb:=proc(f,x,x0,K1,K2) local L,y,i: y:=x0: for i from 1 to K1 do y:=subs(x=y,f): od: L:=[y]: for i from K1+1 to K2 do y:=subs(x=y,f): L:=[op(L),y]: od: evalf(L,10): end: #SS(f,x): Given a function f of x outputs all the fixed points. Try: #SS(5/2*x*(1-x),x); SS:=proc(f,x): [solve(x=f,x)]: end: #SSS(f,x): inputs a function f of x, outputs all the STABLE steady-states. Try: #SSS(5/2*x*(1-x),x); SSS:=proc(f,x) local L,f1,LS,i: L:=SS(f,x): f1:=diff(f,x): LS:=[]: for i from 1 to nops(L) do if abs(evalf(subs(x=L[i],f1)))<1 then LS:=[op(LS),L[i]]: fi: od: LS: end: #Feig(): The list (from wikipedia) of length 8, of the period-doubling cut-offs #of the logistic equation x(n+1)=r*x(n)*(1-x(n)). Calling this list L, i.e. #L:=Feig(): The r below L[1] every orbit tends to ONE ultimate value #for r between L[1] and L[2], in the long run it alternates between 2 ultimate values #between L[2] and L[3] , in the long run it alternates between 4=2^2 ultimate values #Between L[i] and L[i+1], in the long run it alternates between 2^i ultimate values. Try: #Feig(): Feig:=proc() [3,3.4494897,3.5440903,3.5644073,3.5687594,3.5696916,3.5698913,3.5699340]: end: #Digits:=100: #ChaosStory(): This illustarates the period-doubling phenomenon. Try: #ChaosStory() ChaosStory:=proc() local L,orb1,orb2,r,x,f1,r1,i: L:=Feig(): print(`We will investigate the ultimate behaviour of orbits of the discete logistic equation`): print(``): print( ` x(n+1)=r*x(n)*(1-x(n)) `): print(``): print(`as the reproduction paramter r goes up from `): print(``): print(`We will look at what happens after iterating it 10000 times `): print(``): print(`We will have the starting points 0.3 and 0.7 and realize that in the long run it does not matter where you start`): print(``): print(`There are two steady-states, x=0 and x=(r-1)/r `): print(``): print(`The derivative of f(x)=r*x*(1-x) `): print(``): f1:=expand(diff(r*x*(1-x),x)): print(`is `): print(`Plugging in x=0 we get that f'(0) is`): print(subs(x=0,f1)): print(`This is less than 1 in absolute value if r<1, so if r<1 is a stable state `): print(``): print(`Lets try it out for r=0.9 `): print(``): orb1:=Orb(0.9*x*(1-x),x,0.3,10000,10010): orb2:=Orb(0.9*x*(1-x),x,0.7,10000,10010): print(`If you start at 0.3 you get`): lprint(orb1): print(`If you start at 0.7 you get`): lprint(orb2): print(`So 0 indeed the only stable steady-state `): print(`For r>1 x=0 is no longer stable, but for 11 so x=0 is not a steady steady-state `): print(``): print(`So the only stable state for r between 1 and 3 is x=(r-1)/r `): print(``): print(`Lets try it our for first r=2 `): print(``): orb1:=Orb(2*x*(1-x),x,0.3,10000,10400): orb1:=evalf(orb1,10): orb2:=Orb(2*x*(1-x),x,0.7,10000,10400): orb2:=evalf(orb2,10): print(``): print(`If you start at 0.3 you get`): print(``): lprint(orb1): print(``): print(`If you start at 0.7 you get`): print(``): lprint(orb2): print(``): print(`so in the long-run it converges to an ultimate orbit of period`, nops(convert(orb1,set))): for i from 1 to nops(L)-1 do r1:=(L[i]+L[i+1])/2: print(`Let's see how the orbits behave when r is between`, L[i], `and `, L[i+1] ): print(``): print(`Let's take r to be half-way between them, so take r= `, evalf(r1,10)): orb1:=Orb(r1*x*(1-x),x,0.3,10000,10400): orb1:=evalf(orb1,10): orb2:=Orb(r1*x*(1-x),x,0.7,10000,10400): orb2:=evalf(orb2,10): print(``): print(`If you start at 0.3 you get`): print(``): lprint(orb1): print(``): print(`If you start at 0.7 you get`): print(``): lprint(orb2): print(``): print(`so in the long-run it converges to an ultimate orbit of period`, nops(convert(orb1,set))): od: end: #Orbk(k,z,f,INI,K1,K2): Given a positive integer k, a letter (symbol), z, an expression f of z[1], ..., z[k] (representing a multi-variable function of the variables z[1],...,z[k] #a vector INI representing the initial values [x[1],..., x[k]], and (in applications) positive integres K1 and K2, outputs the #values of the sequence starting at n=K1 and ending at n=K2. of the sequence satisfying the difference equation ##x[n+k]=f(x[n+k-1],x[n+k-2],..., x[n]): #This is a generalization to higher-order difference equation of procedure Orb(f,x,x0,K1,K2). For example #Orbk(1,z,5/2*z[1]*(1-z[1]),[0.5],1000,1010); should be the same as #Orb(5/2*z[1]*(1-z[1]),z[1],[0,5],1000,1010); #Try: #Orbk(2,z,(5/4)*z[1]-(3/8)*z[2],[1,2],1000,1010); Orbk:=proc(k,z,f,INI,K1,K2) local L, i,newguy: L:=INI: #We start out with the list of initial values if not (type(k, integer) and type(z,symbol) and type(INI,list) and nops(INI)=k and type(K1,integer) and type(K2,integer) and K1>=1 and K2>=K1) then #checking that the input is OK print(`bad input`): RETURN(FAIL): fi: while nops(L)=0 and K1<=K2) then print(`bad input`): RETURN(FAIL): fi: x1:=x0: for i from 0 to K1-1 do x1:=[seq(subs( {seq(x[i2]=x1[i2],i2=1..nops(x))},F[i1]),i1=1..nops(F))]: od: L:=[x1]: for i from K1 to K2-1 do x1:=[seq(subs( {seq(x[i2]=x1[i2],i2=1..nops(x))},F[i1]),i1=1..nops(F))]: L:=[op(L),expand(x1)]: #we append it to the list od: evalf(L,10): #that's the output end: #RecToT(k,z,f,INI): Given a k-th order recurrence a(n+k)=f(a(n+k-1), ..., a(n)) where f is written in terms of z[1], ..., z[k], and initial conditions #INI=[a(1),..., a(k)], outputs the transformation from F from R^k to R^k such that #if x(n)=[a(n+k-1), ..., a(n)] then x(n+1)= F(x(n)). Try: #RecToT(2,z,(z[1]+z[2])/(z[1]+3*z[2]),[1,2]); RecToT:=proc(k,z,f,INI) local i: [f,seq(z[i],i=1..k-1)],[seq(z[i],i=1..k)],[seq(INI[k+1-i],i=1..k)]: end: #RecToTs(k,z,f): Given a k-th order recurrence a(n+k)=f(a(n+k-1), ..., a(n)) where f is written in terms of z[1], ..., z[k] #if x(n)=[a(n+k-1), ..., a(n)] then x(n+1)= F(x(n)). Try: #RecToTs(2,z,(z[1]+z[2])/(z[1]+3*z[2])); RecToTs:=proc(k,z,f) local i: [f,seq(z[i],i=1..k-1)],[seq(z[i],i=1..k)]: end: #RT(var,K): A random rational transformation of numerator and denominator degrees 1 from R^k to R^k (where k=nops(var), with postiive integer coefficients from 1 to K The inpus are a list of variables x and a posi. integer K #is for generating examples #Try: #RT([x,y],10); RT:=proc(x,K) local ra,i,i1: if not (type(x,list) and type(K, integer) and K>0) then print(`bad input`): RETURM(FAIL): fi: ra:=rand(1..K): #random integer from -K to K [seq((ra()+add(ra()*x[i1],i1=1..nops(x)))/(ra()+add(ra()*x[i1],i1=1..nops(x))), i=1..nops(x))]: end: #RR(var,K): A random rational function with numerator and denominator of degree 1 in the list of variables var, where the ocefficinets are random from 1 to K Try: #RR([x,y,z],10); #is for generating examples #Try: #RR([x,y],10); RR:=proc(x,K) local ra,i,i1: if not (type(x,list) and type(K, integer) and K>0) then print(`bad input`): RETURM(FAIL): fi: ra:=rand(1..K): #random integer from -K to K (ra()+add(ra()*x[i1],i1=1..nops(x)))/(ra()+add(ra()*x[i1],i1=1..nops(x))): end: #SSg(T,x): Given a transormation T in the list of variables, finds all the steady-states (in floating point). Try: #SSg([(1+x)/(3+y),(3+x)/(1+y)],[x,y]); SSg:=proc(T,x) local eq,var,i: if nops(T)<>nops(x) then RETURN(FAIL): fi: eq:={seq(T[i]=x[i],i=1..nops(T))}: var:=[solve(eq,{op(x)})]: evalf([seq(subs(var[i],x),i=1..nops(var))],10): end: #SSr(f,x): Given the underlying function of a recurrence a(n+k)=f(a(n+k-1), ..., a(n)) where f is expressed in terms of x[1],..., x[k] (and k:=nops(k) #Finds all the steady-states. Try: #SSr((3+x)/(1+y),[x,y]); SSr:=proc(f,x) local z,i,sol: sol:=[solve({normal(z-subs({seq(x[i]=z,i=1..nops(x))},f))},{z})]: [seq(subs(sol[i],z),i=1..nops(sol))]: end: #IsDisStable(M): inputs a numeric matrix M (given as a list of lists M) and decides whether all ite eigenvalues have absolute value less than 1 #IsDisStable(Matrix([[1,-1],[-1,1]])); IsDisStable:=proc(M) local Ei1,i: Ei1:=Eigenvalues(evalf(Matrix(M))): evalb(max(seq(abs(Ei1[i]),i=1..nops(M)))<1): end: #JAC(F,x): The Jacobian Matrix (given as a list of lists) of the transformation F in the list of variables x. Try: #JAC([x+y,x^2+y^2],[x,y]): JAC:=proc(F,x) local i,j: if not (type(F,list) and type(x,list) and nops(F)=nops(x)) then print(`Bad input`): RETURN(FAIL): fi: normal([seq([seq(diff(F[i],x[j]),j=1..nops(x))],i=1..nops(F))]): end: #SSSg(F,x): Given a transformation F in the list of variables finds all the STABLE fixed point of the transformation x->F(x), i.e. the set of solutions of #the system {x[1]=F[1], ..., x[k]=F[k]}} that are stable. Try: #SSSg([5/2*x*(1-x)],[x]]); #SSSg([(1+x+y)/(2+3*x+y), (3+x+2*y)/(5+x+3*y)],[x,y]]); SSSg:=proc(F,x) local i, Fi,St,J,J0,pt: if not (type(F,list) and type(x,list) and nops(F)=nops(x)) then print(`bad input`): RETURN(FAIL): fi: Fi:=evalf(SSg(F,x)): #Fi is the set of fixed points in floating-point St:={}: #St is the set of stable fixed points, that starts out empty J:=JAC(F,x): #The general Jacobian in terms of the list of variables x for pt in Fi do #we examine each fixed point, one at a time J0:=subs({seq(x[i]=pt[i],i=1..nops(x))},J): #J0 is the NUMETRICAL Jacobian matrix evaluated at the examined fixed point if IsDisStable(J0) then St:=St union{pt}: #if it is stable we include it fi: od: St: #The output is the set of all the successful fixed points that happened to be stable end: #NicholsonBailey(L,a,c): The discrete-time, double-species dynamical model of Nicholson and Bailey (1935), given by Eqs. (21a)(21b) in Edelstein-Keshet section 3.2 (p. 81) #where the variables are N (hosts) and parasites (P) and the parameters are L (called Lambda there), a, and c #Try: #NicholsonBailey(L,a,c,N,P); #NicholsonBailey(2,0.068,1,N,P); NicholsonBailey:=proc(L,a,c,N,P) [L*N*exp(-a*P),c*N*(1-exp(-a*P))]: end: #NicholsonBaileyM(a,r,K,N,B): The discrete-time, double-species dynamical model of the MODIFIED Nicholson and Bailey model (1935), given by Eqs. (28a)(28b) in Edelstein-Keshet section 3.4 (p. 84) #where the variables are N (hosts) and parasites (P) and the parameters are r and K #Try: #NicholsonBaileyM(r,a,K,N,P); #NicholsonBaileyM(0.5,0.1,14,N,P); NicholsonBaileyM:=proc(r,a,K,N,P) [N*exp(r*(1-N/K)-a*P),N*(1-exp(-a*P))]: end: #HW3(u,v,w): The Hardy-Weinberg unerlying transformation witu (u,v,w), Eqs. (53a,53b,53c) in Edelestein-Keshet Ch. 3 HW3:=proc(u,v,w): [u^2+u*v+(1/4)*v^2, u*v+2*u*w+1/2*v^2+v*w, 1/4*v^2+v*w+w^2]:end: #HW(u,v): The Hardy-Weinberg unerlying transformation witu (u,v,w), Eqs. (53a,53b,53c) in Edelestein-Keshet Ch. 3 using the fact that u+v+w=1 HW:=proc(u,v): expand([u^2+u*v+(1/4)*v^2 , u*v+2*u*(1-u-v)+1/2*v^2+v*(1-u-v)]),[u,v]:end: #HW3g(u,v,w,M): The Hardy-Weinberg unerlying transformation with (u,v,w), GENERALIZED Eqs. with the 3 by 3 matrix M (53a,53b,53c) in Edelestein-Keshet Ch. 3 #Based on Anne Somalwar's solution of the bonus problem from hw15, see the end of #from https://sites.math.rutgers.edu/~zeilberg/Bio21/HW15posted/hw15AnneSomalwar.pdf HW3g:=proc(u,v,w,M) local tot,LI: LI:=[ M[1][1]*u^2+ (M[1][2]+M[2][1])/2*u*v+M[2][2]*(1/4)*v^2, (M[1][2]+M[2][1])/2*u*v+(M[1][3]+M[3][1])*u*w+M[2][2]/2*v^2+ (M[2][3]+M[3][2])/2*v*w, M[2][2]*1/4*v^2+ (M[2][3]+M[3][2])/2* v*w + M[3][3]*w^2]: tot:=LI[1]+LI[2]+LI[3]: [LI[1]/tot,LI[2]/tot,LI[3]/tot]: end: #HWg(u,v,M): The Generalized Hardy-Weinberg unerlying transformation with (u,v), M is the survival matrix. Based on Ann Somalwar's HW3g(u,v,w) (only retain the first two components and replace w by 1-u-v) HWg:=proc(u,v,M) local LI,w: LI:=HW3g(u,v,w,M): normal(subs(w=1-u-v,[LI[1],LI[2]])): end: