> #OK to post #HW20, Timothy Nasralla, 11/15/21 > read "/Users/tan88/OneDrive - Rutgers University/DMB.txt" First Written: Nov. 2021 This is DMB.txt, A Maple package to explore Dynamical models in Biology (both discrete and continuous) accompanying the class Dynamical Models in Biology, Rutgers University. Taught by Dr. Z. (Doron Zeilbeger) The most current version is available on WWW at: http://sites.math.rutgers.edu/~zeilberg/tokhniot/DMB.txt . Please report all bugs to: DoronZeil at gmail dot com . For general help, and a list of the MAIN functions, type "Help();". For specific help type "Help(procedure_name);" ------------------------------ For a list of the supporting functions type: Help1(); For help with any of them type: Help(ProcedureName); ------------------------------ For a list of the functions that give examples of Discrete-time dynamical systems (some famous), type: HelpDDM(); For help with any of them type: Help(ProcedureName); ------------------------------ For a list of the functions continuous-time dynamical systems (some famous) type: HelpCDM(); For help with any of them type: Help(ProcedureName); ------------------------------ > #Q1: Find the equilibrium points, stable equilibrium points, time series plot, > and phase diagram of the SIRS model with varying beta values, nu values, and > gamma values. > #Nu = 2, gamma = 5 G := SIRS(s, i, 0.6/1000,5, 2, 1000) > EquP(G,[s,i]) {[1000., 0.], [3333.333333, -1666.666667]} > SEquP(G,[s,i]) {[1000., 0.]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G, [s,i], [800,200], 0.01, 10) > #As shown in the Phase Diagram and Time Series, the only equilibrium point for > this population would be 1000 susceptible and 0 infected. > G := SIRS(s, i, 1.8/1000,5, 2, 1000) > EquP(G,[s,i]) {[1000., 0.], [1111.111111, -79.36507937]} > SEquP(G,[s,i]) {[1000., 0.]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #Similar to the first beta value, this phase diagram and time series reflect > that the only stable equilibrium is [1000,0]. It is also notable that it took > a longer time to reach said equilibrium. > G := SIRS(s, i, 7.8/1000,5, 2, 1000) > EquP(G,[s,i]) {[256.4102564, 531.1355311], [1000., 0.]} > SEquP(G,[s,i]) {[256.4102564, 531.1355311]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #As shown in this SIRS model, the stable equilibrium point (unlike the first > two) is [256.4102564, 531.1355311] with more infected people, less susceptible > people, and many people being removed from the population as a whole. > #Nu = 3, gamma = 6 G := SIRS(s, i, 0.9/1000,6, 3, 1000) > EquP(G,[s,i]) {[1000., 0.], [3333.333333, -1555.555556]} > SEquP(G,[s,i]) {[1000., 0.]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 2); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #As shown in this SIRS model, the only stable equilibrium point is [1000, 0]. > G := SIRS(s, i, 2.7/1000,6, 3, 1000) > EquP(G,[s,i]) {[1000., 0.], [1111.111111, -74.07407407]} > SEquP(G,[s,i]) {[1000., 0.]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #Just like in the first set of gamma and nu values with the 2nd beta value, > this SIRS model has the stable equilibrium [1000,0] and took longer to reach > this equilibrium than it did with the first beta value. > G := SIRS(s, i, 11.7/1000,6, 3, 1000) > EquP(G,[s,i]) {[256.4102564, 495.7264957], [1000., 0.]} > SEquP(G,[s,i]) {[256.4102564, 495.7264957]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #The equilibrium point for this system is [256.4102564, 495.7264957] which is > reflected in the phase diagram and time series. > #Nu = 4, gamma = 1 G := SIRS(s, i, 1.2/1000, 1, 4, 1000) > EquP(G,[s,i]) {[1000., 0.], [3333.333333, -466.6666667]} > SEquP(G,[s,i]) {[1000., 0.]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #As shown in this SIRS model, the only stable equilibrium point is [1000, 0] > G := SIRS(s,i, 3.6/1000, 1, 4, 1000) > EquP(G,[s,i]) {[1000., 0.], [1111.111111, -22.22222222]} > SEquP(G,[s,i]) {[1000., 0.]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #As seen in the above phase diagrams and phase planes, the only stable > equilibrium point is [1000, 0]. > G := SIRS(s, i, 15.6/1000, 1, 4, 1000) > EquP(G,[s,i]) {[256.4102564, 148.7179487], [1000., 0.]} > SEquP(G,[s,i]) {[256.4102564, 148.7179487]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #As shown in the above phase diagram, starting at [800, 200] leads to a stable > equilibrium at about [256, 148] as reflected in the SEquP function. > #Nu = 7, gamma = 10 G := SIRS(s, i, 2.1/1000, 10, 7, 1000) > EquP(G,[s,i]) {[1000., 0.], [3333.333333, -1372.549020]} > SEquP(G,[s,i]) {[1000., 0.]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #As shown in the above SIRS model, the only stable equilibrium point is [1000, > 0] > G := SIRS(s, i, 6.3/1000, 10, 7, 1000) > EquP(G,[s,i]) {[1000., 0.], [1111.111111, -65.35947712]} > SEquP(G,[s,i]) {[1000., 0.]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #Similar to the first beta value, the stable equilibrium point is [1000, 0] > G := SIRS(s, i, 27.3/1000, 10, 7, 1000) > EquP(G,[s,i]) {[256.4102564, 437.4057315], [1000., 0.]} > SEquP(G,[s,i]) {[256.4102564, 437.4057315]} > TimeSeries(G,[s,i],[800,200],0.01, 10, 1); TimeSeries(G,[s,i],[800,200],0.01, > 10, 2) > PhaseDiag(G,[s,i],[800,200], 0.01, 10) > #Unlike the first two betas, this stable equilibrium point is [254, 437] as > shown in the phase diagram. > #Q2: Use the RandNice function to generate a transformation from (x,y) --> > F(x,y). Find the equilibrium points, the stable equilibrium points, then > compare to the TimeSeries and PhaseDiagram. > F := RandNice([x,y],3) > EquP(F,[x,y]) / [4 1] [8 -1]\ { [5, -7], [-, -], [-, --] } \ [7 7] [5 5 ]/ > SEquP(F,[x,y]) {[5., -7.]} > TimeSeries(F,[x,y],[5.01,-6.9],0.01, 10, 1); > TimeSeries(F,[x,y],[5.01,-6.9],0.01, 10, 2) > PhaseDiag(F,[x,y],[5.01,-6.99],0.01, 10) > #The horizontal asymptote is shown at [1,0] which reflects the stable > equilibrium point found. > F := RandNice([x,y],3) > EquP(F,[x,y]) / [1 4]\ { [-5, 4], [-1, 4], [3, 0], [-, -] } \ [7 7]/ > SEquP(F,[x,y]) {[0.1428571429, 0.5714285714]} > TimeSeries(F,[x,y],[0.15,0.58],0.01, 10, 1); > TimeSeries(F,[x,y],[0.15,0.58],0.01, 10, 2) > PhaseDiag(F,[x,y],[0.15,0.58],0.01, 10) > #As expected, the PhaseDiag function shows the given point returning to the > found equilibrium. > F := RandNice([x,y],3) > EquP(F,[x,y]) /[ 1] [ 3] [-7 8] [7 -3]\ { [0, -], [0, -], [--, -], [-, --] } \[ 3] [ 2] [3 3] [5 5 ]/ > SEquP(F,[x,y]) {[-2.333333333, 2.666666667], [1.400000000, -0.6000000000]} > TimeSeries(F,[x,y],[-2.32,2.67],0.01, 10, 1); > TimeSeries(F,[x,y],[-2.32,2.67],0.01, 10, 2) > PhaseDiag(F,[x,y],[-2.32,2.67],0.01, 10) > TimeSeries(F,[x,y],[1.41,-0.59],0.01, 10, 1); > TimeSeries(F,[x,y],[1.41,-0.59],0.01, 10, 2) > PhaseDiag(F,[x,y],[1.41,-0.59],0.01, 10) > #Both the first and second equilibrium points are approached over time after > starting from a point 0.01 units away from itself. > #Q3: Use the Orbk function to numerically find any stable equilibrium points > of the given system, then use ToSys and SEquP to find them. > #x(n) = (3 + x(n-2) + x(n-3) + x(n-4)) / (1 + x(n-1) + x(n-3)) --> z = (3 + > z[2] + z[3] + z[4]) / (1 + z[1] + z[3]) [z[1] + z[2], z[1]], [z[1], z[2]] [3 + z[2] + z[3] + z[4] ] [----------------------, z[1], z[2], z[3]], [z[1], z[2], z[3], z[4]] [ 1 + z[1] + z[3] ] > Orbk(4,z,(3 + z[2] + z[3] + z[4]) / (1 + z[1] + z[3]), [10., 0., 1., 4.], > 1000, 1010) [1.887962406, 1.760841066, 1.887962406, 1.760841066, 1.887962406, 1.760841066, 1.887962406, 1.760841066, 1.887962406, 1.760841066, 1.887962406] > Orbk(4,z,(3 + z[2] + z[3] + z[4]) / (1 + z[1] + z[3]), [2., 4., 6., 8.], 1000, > 1010) [3.218058756, 1.143841858, 3.218058756, 1.143841858, 3.218058756, 1.143841858, 3.218058756, 1.143841858, 3.218058756, 1.143841858, 3.218058756] > Orbk(4,z,(3 + z[2] + z[3] + z[4]) / (1 + z[1] + z[3]), [1., 1., 1., 1.], 1000, > 1010) [1.567364551, 2.139552296, 1.567364551, 2.139552296, 1.567364551, 2.139552296, 1.567364551, 2.139552296, 1.567364551, 2.139552296, 1.567364551] > Orbk(4,z,(3 + z[2] + z[3] + z[4]) / (1 + z[1] + z[3]), [1.75, 1.75, 1.75, > 1.75], 1000, 1010) [1.803873399, 1.842154845, 1.803873399, 1.842154845, 1.803873399, 1.842154845, 1.803873399, 1.842154845, 1.803873399, 1.842154845, 1.803873399] > #Seemingly there are no stable fixed points for this function. > Q := ToSys(4,z,(3 + z[2] + z[3] + z[4]) / (1 + z[1] + z[3])) > evalf(SFP(Q)); {[1.822875656, 1.822875656, 1.822875656, 1.822875656]} > SEquP(Q,z) {} > Orbk(4,z,(3 + z[2] + z[3] + z[4]) / (1 + z[1] + z[3]), [1.83, 1.83, 1.83, > 1.83], 1000, 1010) [1.824704696, 1.821049140, 1.824704696, 1.821049140, 1.824704696, 1.821049140, 1.824704696, 1.821049140, 1.824704696, 1.821049140, 1.824704696] > #Seemingly the Orbk function of the supposed stable fixed point doesn't > alternate shows that the long term future of the transformation from the given > stable fixed point by SFP isn't actually stable at that point. It is also > worth noting that the SEquP returns the empty set.