#DMB12.txt: Maple Code for Dr. Z.'s Dynamical Models in Biology class, Lecture 12, Oct. 10, 2025

with(linalg):

print(`For a list of the Main procedures type: Help12(); for help with a specific procedure type: Help12(ProcedureName); for example Help12(Feig);`):



Help12a:=proc()
if nargs=0 then
print(`The other procedures are: `):
else
Help12(args):
fi:
end:


Help12:=proc()
if nargs=0 then
print(`The available procedures are: ChaosStory, Feig, Orb, ORB, Orbk, RecToT, SS, SSg, SSS, 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]=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]=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]=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);`):

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 1<r<3 x=(r-1)/r is `):

print(`Indeed plugging-in x=(r-1)/r in f'(x)=`, f1):
print(`You get `):
print(normal(subs(x=(r-1)/r,f1))):
print(`This is less than 1 in absolute value for r from 1 to 3, so x=`, (r-1)/r, `is a stable steady state `):

print(`On the other hand plugging-in x=0 in f'(x)=`, f1):
print(`You get `):
print(normal(subs(x=0,f1))):
print(`Which is larger than 1 for r>1 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)<K2 do
 newguy:=subs({seq(z[i]=L[-i],i=1..k)},f):   #Using what we know about the value yesterday, the day before yesterday, ... up to k days before yesterday we find the value of the sequence today
  L:=[op(L),newguy] :  #we append the new value to the running list of values of our sequence
  od:

[op(K1..K2,L)]:

end:


#ORB(F,x,x0,K1,K2): Inputs a transformation F in the list of variables x with initial point pt, outputs the trajectory
#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.
#For the full trajectory (from n=0 to n=K2), use K1=0. Try:
#ORB([5/2*x*(1-x)],[x], [0.5], 1000,1010);
#ORB([(1+x+y)/(2+x+y),(1+x)/(1+y)],[x,y], [2.,3.], 1000,1010);
ORB:=proc(F,x,x0,K1,K2) local x1,i,L,i1,i2:

if not (type(F,list) and type(x,list) and type(x0,list) and nops(F)=nops(x) and  nops(x)=nops(x0) and type(K1, integer) and type(K2,integer) and K1>=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:


#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:





#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: