CoCalc -- Capstone Project--Epidemiology .ipynb

⁸⁴ views
ubuntu2004

Kernel: SageMath 9.2

In [1]:

#weeks 7&8 
#simulate the model for a range of parameter values and initial conditions 
var('S','A','I','R')
beta=0.09
alpha=0.11
m1=0.06
m2=0.12
m3=0.01
d=0.05
d2=0.03
Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
Aprime=(alpha*S*A)-(m1*A)
Iprime=(beta*S*I)-(m2*I)-(d*I)
Rprime=(m1*A)+(m2*I)-(m3*R)
t=srange(0,100,0.1)
sol=desolve_odeint([Sprime,Aprime,Iprime,Rprime],times=t,ics=[1000,100,150,400],dvars=[S,A,I,R])
p=list_plot(list(zip(t,sol[:,0])),color="purple",legend_label="Susceptable",plotjoined=true)+list_plot(list(zip(t,sol[:,1])),color="green",legend_label="Asymptomatic",plotjoined=true)+list_plot(list(zip(t,sol[:,2])),color="red",legend_label="infected",plotjoined=true)+list_plot(list(zip(t,sol[:,3])),color="blue",legend_label="Recovered",plotjoined=true,axes_labels=["time","population"])
show(p)

Out[1]:

In [2]:

@interact
def icChanger(icS=(0,1000,50),icA=(0,1000,50),icI=(0,1000,50),icR=(0,1000,50)):
    var('S','A','I','R')
    beta = 0.09
    alpha = 0.11
    m1 = 0.06
    m2 = 0.12
    m3 = 0.01
    d = 0.05
    d2=0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    t = srange(0,100,0.1)
    sol = desolve_odeint([Sprime,Aprime,Iprime,Rprime], times = t, ics =
    [icS,icA,icI,icR],dvars=[S,A,I,R])
    p = list_plot(list(zip(t,sol[:,0])),color = "purple",legend_label =
    "Susceptible", plotjoined=true) + list_plot(list(zip(t,sol[:,1])),color =
    "green",legend_label = "Asymptomatic", plotjoined=true) + list_plot(list(zip(t,sol[:,2])),color = "red",legend_label = "Infected",plotjoined=true) + list_plot(list(zip(t,sol[:,3])),color = "blue",legend_label ="Recovered", plotjoined=true, axes_labels = ["t","Pop"])
    show(p)

Out[2]:

Interactive function <function icChanger at 0x7fab240e6a60> with 4 widgets
  icS: IntSlider(value=500, descrip…

In [3]:

@interact
def BetaChanger(beta=(0,1,0.1)):
    var('S','A','I','R')
    alpha = 0.11
    m1 = 0.06
    m2 = 0.12
    m3 = 0.01
    d = 0.05
    d2=0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    t = srange(0,100,0.1)
    sol = desolve_odeint([Sprime,Aprime,Iprime,Rprime], times = t, ics =
    [1000,100,150,400],dvars=[S,A,I,R])
    p = list_plot(list(zip(t,sol[:,0])),color = "purple",legend_label =
    "Susceptible", plotjoined=true) + list_plot(list(zip(t,sol[:,1])),color =
    "green",legend_label = "Asymptomatic", plotjoined=true) + list_plot(list(zip(t,sol[:,2])),color = "red",legend_label = "Infected",plotjoined=true) + list_plot(list(zip(t,sol[:,3])),color = "blue",legend_label ="Recovered", plotjoined=true, axes_labels = ["t","Pop"])
    show(p)

Out[3]:

Interactive function <function BetaChanger at 0x7fa8ff50e9d0> with 1 widget
  beta: FloatSlider(value=0.0, des…

In [4]:

@interact
def AlphaChanger(alpha=(0,1,0.01)):
    var('S','A','I','R')
    beta = 0.09
    m1 = 0.06
    m2 = 0.12
    m3 = 0.01
    d = 0.05
    d2=0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    t = srange(0,100,0.1)
    sol = desolve_odeint([Sprime,Aprime,Iprime,Rprime], times = t, ics =
    [1000,100,150,400],dvars=[S,A,I,R])
    p = list_plot(list(zip(t,sol[:,0])),color = "purple",legend_label =
    "Susceptible", plotjoined=true) + list_plot(list(zip(t,sol[:,1])),color =
    "green",legend_label = "Asymptomatic", plotjoined=true) + list_plot(list(zip(t,sol[:,2])),color = "red",legend_label = "Infected",plotjoined=true) + list_plot(list(zip(t,sol[:,3])),color = "blue",legend_label ="Recovered", plotjoined=true, axes_labels = ["t","Pop"])
    show(p)

Out[4]:

Interactive function <function AlphaChanger at 0x7fa8ff50e820> with 1 widget
  alpha: FloatSlider(value=0.0, d…

In [5]:

@interact
def m1Changer(m1=(0,1,0.01)):
    var('S','A','I','R')
    alpha=0.11
    beta = 0.09
    m2 = 0.12
    m3 = 0.01
    d = 0.05
    d2=0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    t = srange(0,100,0.1)
    sol = desolve_odeint([Sprime,Aprime,Iprime,Rprime], times = t, ics =
    [1000,100,150,400],dvars=[S,A,I,R])
    p = list_plot(list(zip(t,sol[:,0])),color = "purple",legend_label =
    "Susceptible", plotjoined=true) + list_plot(list(zip(t,sol[:,1])),color =
    "green",legend_label = "Asymptomatic", plotjoined=true) + list_plot(list(zip(t,sol[:,2])),color = "red",legend_label = "Infected",plotjoined=true) + list_plot(list(zip(t,sol[:,3])),color = "blue",legend_label ="Recovered", plotjoined=true, axes_labels = ["t","Pop"])
    show(p)

Out[5]:

Interactive function <function m1Changer at 0x7fa8fa9538b0> with 1 widget
  m1: FloatSlider(value=0.0, descrip…

In [6]:

@interact
def m2Changer(m2=(0,1,0.01)):
    var('S','A','I','R')
    alpha=0.11
    beta = 0.09
    m1=0.06
    m3 = 0.01
    d = 0.05
    d2=0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    t = srange(0,100,0.1)
    sol = desolve_odeint([Sprime,Aprime,Iprime,Rprime], times = t, ics =
    [1000,100,150,400],dvars=[S,A,I,R])
    p = list_plot(list(zip(t,sol[:,0])),color = "purple",legend_label =
    "Susceptible", plotjoined=true) + list_plot(list(zip(t,sol[:,1])),color =
    "green",legend_label = "Asymptomatic", plotjoined=true) + list_plot(list(zip(t,sol[:,2])),color = "red",legend_label = "Infected",plotjoined=true) + list_plot(list(zip(t,sol[:,3])),color = "blue",legend_label ="Recovered", plotjoined=true, axes_labels = ["t","Pop"])
    show(p)

Out[6]:

Interactive function <function m2Changer at 0x7fa900de71f0> with 1 widget
  m2: FloatSlider(value=0.0, descrip…

In [7]:

@interact
def m3Changer(m3=(0,1,0.01)):
    var('S','A','I','R')
    beta = 0.09
    alpha=0.11
    m1=0.06
    m2 = 0.12
    d = 0.05
    d2=0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    t = srange(0,100,0.1)
    sol = desolve_odeint([Sprime,Aprime,Iprime,Rprime], times = t, ics =
    [1000,100,150,400],dvars=[S,A,I,R])
    p = list_plot(list(zip(t,sol[:,0])),color = "purple",legend_label =
    "Susceptible", plotjoined=true) + list_plot(list(zip(t,sol[:,1])),color =
    "green",legend_label = "Asymptomatic", plotjoined=true) + list_plot(list(zip(t,sol[:,2])),color = "red",legend_label = "Infected",plotjoined=true) + list_plot(list(zip(t,sol[:,3])),color = "blue",legend_label ="Recovered", plotjoined=true, axes_labels = ["t","Pop"])
    show(p)

Out[7]:

Interactive function <function m3Changer at 0x7fa8fa54c8b0> with 1 widget
  m3: FloatSlider(value=0.0, descrip…

In [8]:

@interact
def dChanger(d=(0,1,0.01)):
    var('S','A','I','R')
    beta = 0.09
    alpha=0.11
    m1=0.06
    m2 = 0.12
    m3 = 0.01
    d2=0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    t = srange(0,100,0.1)
    sol = desolve_odeint([Sprime,Aprime,Iprime,Rprime], times = t, ics =
    [1000,100,150,400],dvars=[S,A,I,R])
    p = list_plot(list(zip(t,sol[:,0])),color = "purple",legend_label =
    "Susceptible", plotjoined=true) + list_plot(list(zip(t,sol[:,1])),color =
    "green",legend_label = "Asymptomatic", plotjoined=true) + list_plot(list(zip(t,sol[:,2])),color = "red",legend_label = "Infected",plotjoined=true) + list_plot(list(zip(t,sol[:,3])),color = "blue",legend_label ="Recovered", plotjoined=true, axes_labels = ["t","Pop"])
    show(p)

Out[8]:

Interactive function <function dChanger at 0x7fa8fa54c3a0> with 1 widget
  d: FloatSlider(value=0.0, descripti…

In [9]:

#week 9 finding the EPs of the model. We added a new parameter d2 to represent the natural death rate of the susceptible population
var('S','A','I','R')
beta=0.09
alpha=0.11
m1=0.06
m2=0.12
m3=0.01
d=0.05
d2=0.03
Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
Aprime = (alpha*S*A)-(m1*A)
Iprime = (beta*S*I)-(m2*I)-(d*I)
Rprime = (m1*A)+(m2*I)-(m3*R)
solve([Sprime,Aprime,Iprime,Rprime],[S,A,I,R])

Out[9]:

[[S == 0, A == 0, I == 0, R == 0], [S == (17/9), A == 0, I == (-17/15), R == (-68/5)]]

In [10]:

#week 9 testing the stability of the EP 0,0,0,0)
# we will not be testing the stability of the other EP (17/9,0,-17/5,-68/5) because it involves negative values that do not make sense in the context (population) of our model 
k=jacobian([Sprime,Aprime,Iprime,Rprime],[S,A,I,R])
show(k)

Out[10]:

\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} -0.110000000000000 \, A - 0.0900000000000000 \, I - 0.0300000000000000 & -0.110000000000000 \, S & -0.0900000000000000 \, S & 0.0100000000000000 \\ 0.110000000000000 \, A & 0.110000000000000 \, S - 0.0600000000000000 & 0 & 0 \\ 0.0900000000000000 \, I & 0 & 0.0900000000000000 \, S - 0.170000000000000 & 0 \\ 0 & 0.0600000000000000 & 0.120000000000000 & -0.0100000000000000 \end{array}\right)

In [11]:

#week 9
kmatrix=k(S=0,A=0,I=0,R=0)
show (kmatrix)

Out[11]:

\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} -0.0300000000000000 & 0 & 0 & 0.0100000000000000 \\ 0 & -0.0600000000000000 & 0 & 0 \\ 0 & 0 & -0.170000000000000 & 0 \\ 0 & 0.0600000000000000 & 0.120000000000000 & -0.0100000000000000 \end{array}\right)

In [12]:

#week 9
eig_vals=kmatrix.eigenvalues()
show(eig_vals)

Out[12]:

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{17}{100}, -\frac{3}{100}, -\frac{1}{100}, -\frac{3}{50}\right]

In [13]:

#week 9
#since all the eigenvalues are negative, the EP (0,0,0,0) is a stable node

In [14]:

#week 9
k2matrix=k(S=(17/9),A=0,I=(-17/15),R=(-68/5))
#show (k2matrix)

In [35]:

eig_valss=k2matrix.eigenvalues()
show(eig_valss)

Out[35]:

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{500} \, \sqrt{654} - \frac{27}{500}, \frac{1}{500} \, \sqrt{654} - \frac{27}{500}, \frac{17}{100}, \frac{133}{900}\right]

In [17]:

#week 10  
@interact
def m3Changer(m3Slider=(0,1,0.01)):
    var('S','A','I','R')
    beta = 0.09
    alpha = 0.11
    m1 = 0.06
    m2 = 0.12
    m3 = m3Slider
    d1 = 0.05
    d2 = 0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d1*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    k = jacobian([Sprime,Aprime,Iprime,Rprime], [S,A,I,R])
    kmatrix = k(S=(17/9),A=0,I=(-17/15),R=(-68/5))
    eig_vals = kmatrix.eigenvalues()
    #show(eig_vals) #how to extract only the real part of the eigenvalues? 
    for i in [0,1,2,3]:
        eig_vals[i].n().real()
        
        print(eig_vals[i].n().real())

Out[17]:

Interactive function <function m3Changer at 0x7fa8fa1b4790> with 1 widget
  m3Slider: FloatSlider(value=0.0, d…

In [18]:

#week10 
@interact
def m3Changer(m3=(0,1,0.01)):
    var('S','A','I','R')
    beta = 0.09
    alpha=0.11
    m1=0.06
    m2 = 0.12
    d = 0.05
    d2=0.03
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Aprime = (alpha*S*A)-(m1*A)
    Iprime = (beta*S*I)-(m2*I)-(d*I)
    Rprime = (m1*A)+(m2*I)-(m3*R)
    t = srange(0,100,0.1)
    sol = desolve_odeint([Sprime,Aprime,Iprime,Rprime], times = t, ics =
    [1000,100,150,400],dvars=[S,A,I,R])
    p = list_plot(list(zip(t,sol[:,0])),color = "purple",legend_label =
    "Susceptible", plotjoined=true) + list_plot(list(zip(t,sol[:,1])),color =
    "green",legend_label = "Asymptomatic", plotjoined=true) + list_plot(list(zip(t,sol[:,2])),color = "red",legend_label = "Infected",plotjoined=true) + list_plot(list(zip(t,sol[:,3])),color = "blue",legend_label ="Recovered", plotjoined=true, axes_labels = ["t","Pop"])
    show(p)

Out[18]:

Interactive function <function m3Changer at 0x7fa8fa1b4310> with 1 widget
  m3: FloatSlider(value=0.0, descrip…

In [19]:

#final week
var('S','J','A','I','R')
beta=0.09
alpha=0.11
gamma=0.05
m1=0.06
m2=0.12
m3=0.01
m4=0.2
m5=0.3
d=0.05
d2=0.03
d3=0.01
Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
Jprime=(m4*R)-(gamma*J*I)-(d3*J)
Aprime=(alpha*S*A)-(m1*A)
Iprime=(beta*S*I)+(gamma*J*I)-(m2*I)-(d*I)-(m5*I)
Rprime=(m1*A)+(m2*I)-(m3*R)-(m4*R)+(m5*I)
t=srange(0,100,0.1)
sol=desolve_odeint([Sprime,Jprime,Aprime,Iprime,Rprime],times=t,ics=[1000,500,100,150,400],dvars=[S,J,A,I,R])
p=list_plot(list(zip(t,sol[:,0])),color="purple",legend_label="Susceptable",plotjoined=true)+list_plot(list(zip(t,sol[:,1])),color="green",legend_label="Juveniles",plotjoined=true)+list_plot(list(zip(t,sol[:,2])),color="red",legend_label="Asymptomatic",plotjoined=true)+list_plot(list(zip(t,sol[:,3])),color="blue",legend_label="Infected",plotjoined=true)+ list_plot(list(zip(t,sol[:,4])),color="yellow",legend_label="Recovered",plotjoined=true,axes_labels=["time","population"])
show(p)

Out[19]:

In [20]:

@interact
def gammaChanger(gamma=(0,3,0.01)):
    var('S','J','A','I','R')
    beta=0.09
    alpha=0.11
    m1=0.06
    m2=0.12
    m3=0.01
    m4=0.2
    m5=0.3
    d=0.05
    d2=0.03
    d3=0.01
    Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
    Jprime=(m4*R)33-(gamma*J*I)-(d3*J)
    Aprime=(alpha*S*A)-(m1*A)
    Iprime=(beta*S*I)+(gamma*J*I)-(m2*I)-(d*I)-(m5*I)
    Rprime=(m1*A)+(m2*I)-(m3*R)-(m4*R)+(m5*I)
    t=srange(0,100,0.1)
    sol=desolve_odeint([Sprime,Jprime,Aprime,Iprime,Rprime],times=t,ics=[1000,500,100,150,400],dvars=[S,J,A,I,R])
    p=list_plot(list(zip(t,sol[:,0])),color="purple",legend_label="Susceptable",plotjoined=true)+list_plot(list(zip(t,sol[:,1])),color="green",legend_label="Juveniles",plotjoined=true)+list_plot(list(zip(t,sol[:,2])),color="red",legend_label="Asymptomatic",plotjoined=true)+list_plot(list(zip(t,sol[:,3])),color="blue",legend_label="Infected",plotjoined=true)+ list_plot(list(zip(t,sol[:,4])),color="yellow",legend_label="Recovered",plotjoined=true,axes_labels=["time","population"])
    show(p)

Out[20]:

Interactive function <function gammaChanger at 0x7fa8fa95db80> with 1 widget
  gamma: FloatSlider(value=1.0, d…

In [21]:

#final week
var('S','J','A','I','R')
beta=0.09
alpha=0.11
gamma=0.05
m1=0.06
m2=0.12
m3=0.01
m4=0.2
m5=0.3
d=0.05
d2=0.03
d3=0.01
Sprime = (m3*R)-(beta*S*I)-(alpha*S*A)-(d2*S)
Jprime=(m4*R)-(gamma*J*I)-(d3*J)
Aprime=(alpha*S*A)-(m1*A)
Iprime=(beta*S*I)+(gamma*J*I)-(m2*I)-(d*I)-(m5*I)
Rprime=(m1*A)+(m2*I)-(m3*R)-(m4*R)+(m5*I)
show(solve([Sprime,Jprime,Aprime,Iprime,Rprime],[S,J,A,I,R]))

Out[21]:

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[\left[S = 0, J = 0, A = 0, I = 0, R = 0\right], \left[S = \left(\frac{6}{11}\right), J = \left(-\frac{18}{11}\right), A = \left(-\frac{63}{220}\right), I = 0, R = \left(-\frac{9}{110}\right)\right], \left[S = -\frac{5}{9} \, \sqrt{37} - \frac{28}{9}, J = \sqrt{37} + 15, A = 0, I = \frac{2}{15} \, \sqrt{37} - \frac{17}{15}, R = \frac{4}{15} \, \sqrt{37} - \frac{34}{15}\right], \left[S = \frac{5}{9} \, \sqrt{37} - \frac{28}{9}, J = -\sqrt{37} + 15, A = 0, I = -\frac{2}{15} \, \sqrt{37} - \frac{17}{15}, R = -\frac{4}{15} \, \sqrt{37} - \frac{34}{15}\right], \left[S = \left(\frac{6}{11}\right), J = \left(\frac{463}{55}\right), A = \left(\frac{44611}{60500}\right), I = \left(-\frac{553}{275}\right), R = \left(-\frac{115287}{30250}\right)\right]\right]

In [22]:

#final week
newwk=jacobian([Sprime,Jprime,Aprime,Iprime,Rprime],[S,J,A,I,R])
show(newwk)

Out[22]:

\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrr} -0.110000000000000 \, A - 0.0900000000000000 \, I - 0.0300000000000000 & 0 & -0.110000000000000 \, S & -0.0900000000000000 \, S & 0.0100000000000000 \\ 0 & -0.0500000000000000 \, I - 0.0100000000000000 & 0 & -0.0500000000000000 \, J & 0.200000000000000 \\ 0.110000000000000 \, A & 0 & 0.110000000000000 \, S - 0.0600000000000000 & 0 & 0 \\ 0.0900000000000000 \, I & 0.0500000000000000 \, I & 0 & 0.0500000000000000 \, J + 0.0900000000000000 \, S - 0.470000000000000 & 0 \\ 0 & 0 & 0.0600000000000000 & 0.420000000000000 & -0.210000000000000 \end{array}\right)

In [23]:

#1st EP
kkmatrix=newwk(S=0,J=0,A=0,I=0,R=0)
eig_vals=kkmatrix.eigenvalues()
show(eig_vals)

Out[23]:

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{47}{100}, -\frac{21}{100}, -\frac{3}{100}, -\frac{1}{100}, -\frac{3}{50}\right]

In [24]:

#the EP (0,0,0,0,0) is a stable node

In [25]:

#*1st EP (showing only the real part and approximating)
kkmatrix=newwk(S=0,J=0,A=0,I=0,R=0)
eig_vals=kkmatrix.eigenvalues()
for i in [0,1,2,3]:
        eig_vals[i].n().real()
        
        print(eig_vals[i].n().real())

Out[25]:

-0.470000000000000
-0.210000000000000
-0.0300000000000000
-0.0100000000000000

In [36]:

#2nd EP
kkmatrix=newwk(S=(6/11),J=(-18/11),A=(-63/220),I=(0),R=(-9/110))
eig_vals=kkmatrix.eigenvalues()
show(eig_vals)

Out[36]:

\renewcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, {\left(\frac{3}{800000000} i \, \sqrt{8710571163} - \frac{1786609}{8000000000}\right)}^{\frac{1}{3}} {\left(i \, \sqrt{3} + 1\right)} - \frac{22261 \, {\left(-i \, \sqrt{3} + 1\right)}}{8000000 \, {\left(\frac{3}{800000000} i \, \sqrt{8710571163} - \frac{1786609}{8000000000}\right)}^{\frac{1}{3}}} - \frac{139}{2000}, -\frac{1}{2} \, {\left(\frac{3}{800000000} i \, \sqrt{8710571163} - \frac{1786609}{8000000000}\right)}^{\frac{1}{3}} {\left(-i \, \sqrt{3} + 1\right)} - \frac{22261 \, {\left(i \, \sqrt{3} + 1\right)}}{8000000 \, {\left(\frac{3}{800000000} i \, \sqrt{8710571163} - \frac{1786609}{8000000000}\right)}^{\frac{1}{3}}} - \frac{139}{2000}, {\left(\frac{3}{800000000} i \, \sqrt{8710571163} - \frac{1786609}{8000000000}\right)}^{\frac{1}{3}} + \frac{22261}{4000000 \, {\left(\frac{3}{800000000} i \, \sqrt{8710571163} - \frac{1786609}{8000000000}\right)}^{\frac{1}{3}}} - \frac{139}{2000}, -\frac{553}{1100}, -\frac{1}{100}\right]

In [27]:

#the EP ((6/11),(-18/11),(-63/220),0,(-9/110)) is a saddle point

In [28]:

#*2nd EP (approximating)
kkmatrix=newwk(S=(6/11),J=(-18/11),A=(-63/220),I=(0),R=(-9/110))
eig_vals=kkmatrix.eigenvalues()
for i in [0,1,2,3]:
        eig_vals[i].n()
        
        print(eig_vals[i].n())

Out[28]:

-0.0414216546529095 + 2.77555756156289e-17*I
-0.210442501597286 - 1.38777878078145e-17*I
0.0433641562501955 - 1.38777878078145e-17*I
-0.502727272727273

In [29]:

#3rd EP
kkmatrix=newwk(S=(-(5/9)*sqrt(37)-(28/9)),J=sqrt(37)+15,A=0,I=((2/15)*sqrt(37)-(17/15)),R=((4/15)*sqrt(37)-(34/15)))
eig_vals=kkmatrix.eigenvalues()
show(eig_vals)
#this EP is a saddle point

Out[29]:

ParseError: KaTeX parse error: Too many expansions: infinite loop or need to increase maxExpand setting

In [30]:

#*3rd EP (showing only the real part and approximating)
kkmatrix=newwk(S=(-(5/9)*sqrt(37)-(28/9)),J=sqrt(37)+15,A=0,I=((2/15)*sqrt(37)-(17/15)),R=((4/15)*sqrt(37)-(34/15)))
eig_vals=kkmatrix.eigenvalues()
for i in [0,1,2,3]:
        eig_vals[i].n()
        
        print(eig_vals[i].n())

Out[30]:

-0.235733076179253
0177133099044273
00657076618796292 - 0.0766326752817614*I
00657076618796293 + 0.0766326752817614*I

In [31]:

#4th EP
kkmatrix=newwk(S=((5/9)*sqrt(37)-(28/9)),J=-1*sqrt(37)+15,A=0,I=((-2/15)*sqrt(37)-(17/15)),R=((-4/15)*sqrt(37)-(34/15)))
eig_vals=kkmatrix.eigenvalues()
show(eig_vals)
#this EP is a saddle point

Out[31]:

ParseError: KaTeX parse error: Too many expansions: infinite loop or need to increase maxExpand setting

In [32]:

#*4th EP (only approximating values)
kkmatrix=newwk(S=((5/9)*sqrt(37)-(28/9)),J=-1*sqrt(37)+15,A=0,I=((-2/15)*sqrt(37)-(17/15)),R=((-4/15)*sqrt(37)-(34/15)))
eig_vals=kkmatrix.eigenvalues()
for i in [0,1,2,3]:
        eig_vals[i].n()
        
        print(eig_vals[i].n())

Out[32]:

-0.317833294503382 + 2.52299236050494e-18*I
-0.0156516448336855 + 1.11439536231446e-17*I
0.136939385115158 - 1.24986920702229e-17*I
0.218757121454142 - 1.16825391342659e-18*I

In [33]:

#*5th EP (showing only the real part and approximating)
kkmatrix=newwk(S=(6/11),J=(463/55),A=(44611/60500),I=(-553/275),R=(-115287/30250))
eig_vals=kkmatrix.eigenvalues()
for i in [0,1,2,3]:
        eig_vals[i].n()
        
        print(eig_vals[i].n())

Out[33]:

---------------------------------------------------------------------------
ArithmeticError                           Traceback (most recent call last)
<ipython-input-33-f5684b383291> in <module>
      1 #*5th EP (showing only the real part and approximating)
      2 kkmatrix=newwk(S=(Integer(6)/Integer(11)),J=(Integer(463)/Integer(55)),A=(Integer(44611)/Integer(60500)),I=(-Integer(553)/Integer(275)),R=(-Integer(115287)/Integer(30250)))
----> 3 eig_vals=kkmatrix.eigenvalues()
      4 for i in [Integer(0),Integer(1),Integer(2),Integer(3)]:
      5         eig_vals[i].n()
/ext/sage/sage-9.2/local/lib/python3.8/site-packages/sage/matrix/matrix_symbolic_dense.pyx in sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense.eigenvalues (build/cythonized/sage/matrix/matrix_symbolic_dense.c:2937)()
    181         maxima_evals = self._maxima_(maxima).eigenvalues()._sage_()
    182         if not len(maxima_evals):
--> 183             raise ArithmeticError("could not determine eigenvalues exactly using symbolic matrices; try using a different type of matrix via self.change_ring(), if possible")
    184         return sum([[ev] * int(mult) for ev, mult in zip(*maxima_evals)], [])
    185 
ArithmeticError: could not determine eigenvalues exactly using symbolic matrices; try using a different type of matrix via self.change_ring(), if possible

In [34]:

#5th EP
kkmatrix=newwk(S=(6/11),J=(463/55),A=(44611/60500),I=(-553/275),R=(-115287/30250))
eig_vals=kkmatrix.eigenvalues()
show(eig_vals)
#the eigenvalues do not give helpful information for this set of EP. Therefore, we will neglect this EP.

Out[34]:

---------------------------------------------------------------------------
ArithmeticError                           Traceback (most recent call last)
<ipython-input-34-5626d3d66e4e> in <module>
      1 #5th EP
      2 kkmatrix=newwk(S=(Integer(6)/Integer(11)),J=(Integer(463)/Integer(55)),A=(Integer(44611)/Integer(60500)),I=(-Integer(553)/Integer(275)),R=(-Integer(115287)/Integer(30250)))
----> 3 eig_vals=kkmatrix.eigenvalues()
      4 show(eig_vals)
      5 #the eigenvalues do not give helpful information for this set of EP. Therefore, we will neglect this EP.
/ext/sage/sage-9.2/local/lib/python3.8/site-packages/sage/matrix/matrix_symbolic_dense.pyx in sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense.eigenvalues (build/cythonized/sage/matrix/matrix_symbolic_dense.c:2937)()
    181         maxima_evals = self._maxima_(maxima).eigenvalues()._sage_()
    182         if not len(maxima_evals):
--> 183             raise ArithmeticError("could not determine eigenvalues exactly using symbolic matrices; try using a different type of matrix via self.change_ring(), if possible")
    184         return sum([[ev] * int(mult) for ev, mult in zip(*maxima_evals)], [])
    185 
ArithmeticError: could not determine eigenvalues exactly using symbolic matrices; try using a different type of matrix via self.change_ring(), if possible

In [0]:

In [0]:

Product

Resources

Company