#M13.txt: Maple code for Lecture 13 of Dynamical Modesl in Biology, Fall 2021 (taught by Dr. Z.) Help13:=proc(): print(` RT2(x,y,d,K), Orb2(F,x,y,pt0,K1,K2), FP2(F,x,y), SFP2(F,x,y), PlotOrb2(L), FP2drz(F,x,y), SFP2drz(F,x,y) `):end: with(LinearAlgebra): #RT2(x,y,d,K): A random rational transformation of degree d from R^2 to R^2 with postiive integer coefficients from 1 to K The inputs are variables x and y and #the output is a pair of expressions of (x,y) representing functions. It is for generating examples #Try: #RT2(x,y,2,10); RT2:=proc(x,y,d,K) local ra,i,j,f,g: ra:=rand(1..K): #random integer from -K to K f:=add(add(ra()*x^i*y^j,j=0..d-i),i=0..d)/add(add(ra()*x^i*y^j,j=0..d-i),i=0..d): g:=add(add(ra()*x^i*y^j,j=0..d-i),i=0..d)/add(add(ra()*x^i*y^j,j=0..d-i),i=0..d): [f,g]: end: #Orb2(F,x,y,pt,K1,K2): Inputs a mapping F=[f,g] from R^2 to R^2 where f and g describe functions of x and y, an initial point pt0=[x0,y0] #outputs the orbit starting at discrete time K1 and ending in discrete time K2. Try #F:=RT2(x,y,2,10); #Orb2(F,x,y,[1.1,1.2],1000,1010); Orb2:=proc(F,x,y,pt0,K1,K2) local pt,L,i: pt:=pt0: for i from 1 to K1 do pt:=subs({x=pt[1],y=pt[2]},F): od: L:=[]: for i from K1+1 to K2 do L:=[op(L),pt]: pt:=subs({x=pt[1],y=pt[2]},F): od: L: end: #FP2(F,x,y): The list of fixed points of the transformation [x,y]->F. Try #FP2([x-y,x=y],x,y); FP2:=proc(F,x,y) local L,i: L:=[solve({F[1]=x,F[2]=y},{x,y})]: [seq(subs(L[i],[x,y]),i=1..nops(L))]: end: #SFP2(F,x,y): The list of Stable fixed points of the transformation [x,y]->F. Try #SFP2([(1+x)/(1+y), (1+7*y)/(4+x)],x,y); SFP2:=proc(F,x,y) local L,J,S,J0,i,pt,EV: L:=evalf(FP2(F,x,y)): #F is the list of ALL fixed points of the transformation [x,y]->F using the previous procedure FP2(F,x,y), but since we are interested in numbers we take the floating point version using evalf J:=Matrix(normal([[diff(F[1],x),diff(F[1],y)],[diff(F[2],x),diff(F[2],y)]])): #J is the Jacobian matrix in general (in terms of the variables x and y). Note that J is a SYMBOLIC matrix featuring variables x and y S:=[]: #S is the list of stable fixed points that starts out empty for i from 1 to nops(L) do #we examime it case by case pt:=L[i]: #pt is the current fixed point to be examined J0:=subs({x=pt[1],y=pt[2]},J): #J0 is the NUMERICAL matrix obtained by plugging-in the examined fixed pt EV:=Eigenvalues(J0): # We used Maple's command Eigenavalues to find the eigenvalues of this 2 by 2 matrix if abs(EV[1])<1 and abs(EV[2])<1 then S:=[op(S),pt]: #If both eigenvalues have absolute value less than 1 it means that they are stable, so we append the examined fixed point, pt, to the list of fixed points fi: od: S: #the output is S end: ###added Oct. 17, 20221 with(plots): PlotOrb1:=proc(L) local i,d: d:=textplot([L[1],0,0]): for i from 2 to nops(L) do d:=d,textplot([L[i],0,i-1]): od: display(d): end: PlotOrb2:=proc(L) local i,d: d:=textplot([op(L[1]),0]): for i from 2 to nops(L) do d:=d,textplot([op(L[i]),i-1]): od: display(d): end: ###End added Oct. 17, 20221 ###old stuff #M11.txt: Maple code for Lecture 11 of Dynamical Models in Biology taught by Dr. Z. Help11:=proc(): print(` SFPe(f,x), Orbk(k,z,f,INI,K1,K2) `):end: #SFPe(f,x): The set of fixed points of x->f(x) done exactly (and allowing symbolic parameters), followed by the condition of stability (if it is netween -1 and 1 it is stable) #Try: FPe(k*x*(1-x),x); #VERSION OF Oct. 12, 2021 (avoiding division by 0) SFPe:=proc(f,x) local f1,L,i,M: f1:=normal(diff(f,x)): L:=[solve(numer(f-x),x)]: M:=[]: for i from 1 to nops(L) do if subs(x=L[i],denom(f1) )<>0 then M:=[op(M),[L[i],normal(subs(x=L[i],f1))]]: fi: od: M: end: #Added after class #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]=f(x[n-1],x[n-2],..., x[n-k+1]): #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>0 and K2>K1) then #checking that the input is OK print(`bad input`): RETURN(FAIL): fi: while nops(L)f where f is an expression in x. Try: #FP(2*x*(1-x),x); FP:=proc(f,x) evalf([solve(f=x,x)]): end: #SFP(f,x): The list of stable fixed points of the map x->f where f is an expression in x. Try: #SFP(2*x*(1-x),x); SFP:=proc(f,x) local L,i,f1,pt,Ls: L:=FP(f,x): #The list of fixed points (including complex ones) Ls:=[]: #Ls is the list of stable fixed points, that starts out as the empty list f1:=diff(f,x): #The derivative of the function f w.r.t. x for i from 1 to nops(L) do pt:=L[i]: if abs(subs(x=pt,f1))<1 then Ls:=[op(Ls),pt]: # if pt, is stable we add it to the list of stable points fi: od: Ls: #The last line is the output end: #Comp(f,x): f(f(x)) Comp:=proc(f,x): normal(subs(x=f,f)):end: ##added Oct. 17, 2021 #FP2drz(F,x,y): The list of fixed points of the transformation [x,y]->F. Dr. Z.'s way #FP2([x-y,x+y],x,y); FP2drz:=proc(F,x,y) local eq,i,L,S1: eq:=[numer(F[1]-x),numer(F[2]-y)]: L:=Groebner[Basis](eq,plex(x,y)): S1:=evalf([solve(L[1],y)]): [seq([solve(subs(y=S1[i],L[2]),x),S1[i]],i=1..nops(S1))]: end: #SFP2drz(F,x,y): The list of Stable fixed points of the transformation [x,y]->F. Try #SFP2drz([(1+x)/(1+y), (1+7*y)/(4+x)],x,y); SFP2drz:=proc(F,x,y) local L,J,S,J0,i,pt,EV: L:=FP2drz(F,x,y): #F is the list of ALL fixed points of the transformation [x,y]->F using the previous procedure FP2(F,x,y), but since we are interested in numbers we take the floating point version using evalf J:=Matrix(normal([[diff(F[1],x),diff(F[2],x)],[diff(F[1],y),diff(F[2],y)]])): #J is the Jacobian matrix in general (in terms of the variables x and y). Note that J is a SYMBOLIC matrix featuring variables x and y S:=[]: #S is the list of stable fixed points that starts out empty for i from 1 to nops(L) do #we examime it case by case pt:=L[i]: #pt is the current fixed point to be examined J0:=subs({x=pt[1],y=pt[2]},J): #J0 is the NUMERICAL matrix obtained by plugging-in the examined fixed pt EV:=Eigenvalues(J0): # We used Maple's command Eigenavalues to find the eigenvalues of this 2 by 2 matrix if abs(EV[1])<1 and abs(EV[2])<1 then S:=[op(S),pt]: #If both eigenvalues have absolute value less than 1 it means that they are stable, so we append the examined fixed point, pt, to the list of fixed points fi: od: S: #the output is S end: