#Maple code for Lecture 18 Help18:=proc(): print(` Dis2(F,x,y,pt,h,A), SIRS(s,i,beta,gamma,nu,N) `):end: #SIRS(s,i,beta,gamma,nu,N): The SIRS dynamical model with parameters beta,gamma, nu,N (see section 6.6 of Edelstein-Keshet), s is the number of #Susceptibles, i is the number of infected, (the number of removed is given by N-s-i). N is the total population SIRS:=proc(s,i,beta,gamma,nu,N):[-beta*s*i+gamma*(N-s-i),beta*s*i-nu*i]:end: #Dis2(F,x,y,pt,h,A): The approximate orbit of the Dynamical system approximating the 2D for the autonomous continuous dynamical process #dx/dt=F[1](x(t),y(t)) #dy/dt=F[2](x(t),y(t)) , x(0)=pt[1], y(0)=pt[2] with mesh size h from t=0 to t=A Dis2:=proc(F,x,y,pt,h,A) local L,i: L:=Orb2([x+h*F[1],y+h*F[2]], x,y, pt,0,trunc(A/h)) : L:=[seq([i*h,[L[i][1],L[i][2]]],i=1..nops(L))]: end: #OLD STUFF Help17:=proc(): print(` HW3g(u,v,w,M), HW2g(u,v,M) `):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: #HW2g(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) HW2g:=proc(u,v,M) local LI,w: LI:=HW3g(u,v,w,M): normal(subs(w=1-u-v,[LI[1],LI[2]])): end: #OLD STUFF Help15:=proc(): print(` HW3(u,v,w), HW2(u,v) , Dis1(F,x,x0,h,A), ToSys(k,z,f,INI) `):end: #ToSys(k,z,f,INI): converts the kth order difference equation x(n)=f(x[n-1],x[n-2],...x[n-k]) to a first-order system #x1(n)=F(x1(n-1),x2(n-1), ...,xk(n-1)) #x2(n)=x1(n-1) #... #xk(n)=x[k-1](n-1). It gives the underlying transformation phrased in terms of z[1],...z[k], followed by the initial conditions. Try: #ToSys:=proc(2,z,z[1]+z[2],[1,1]) ToSys:=proc(k,z,f,INI) local i: [f,seq(z[i-1],i=2..k)],INI: 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: #HW2(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 HW2:=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)]):end: #Dis1(F,x,x0,h,A): The approximate orbit of the Dynamical system approximating the 1D for the autonomous continuous dynamical process dy/dt=F(y(t)) , y(0)=y0 with mesh size h from t=0 to t=A Dis1:=proc(F,x,x0,h,A) local L,i: L:=Orb(x+h*F,x,x0,0,trunc(A/h)): L:=[seq([i*h,L[i]],i=1..nops(L))]: end: ##old stuff #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-1 do pt:=subs({x=pt[1],y=pt[2]},F): od: L:=[]: for i from K1 to K2 do L:=[op(L),pt]: pt:=normal(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: