#p is probability of expressed schizophrenic at time t
#q is probability of carrier schizophrenia at time t
var('p q')
(p, q) (p, q) |
#weighted coupling and birth rates
var('nn nc ns cn cc cs sn sc ss norm')
(nn, nc, ns, cn, cc, cs, sn, sc, ss, norm) (nn, nc, ns, cn, cc, cs, sn, sc, ss, norm) |
nn=.934
nc=cn=.86
ns=sn=.77
cc=.68
sc=cs=.59
ss=.50
|
|
#the normalizer
norm(p, q) = nn*(1-p-q)^2 + (nc+cn)*q*(1-p-q) + (ns+sn)*p*(1-p-q) + cc*q^2 + (sc+cs)*p*q + ss*p^2
|
|
#this is the number of people with expressed schizophrenia at time t+1
Pt(p, q) =(.011*nn*(1-p-q)^2 + .07*(nc+cn)*q*(1-p-q) + .13*(ns+sn)*p*(1-p-q) + .16*cc*q^2 + .255*(sc+cs)*p*q + .38*ss*p^2)/norm(p, q)
|
|
#this is the number of people who are carriers of schizophrenia at time t+1
Qt(p, q)=(.027*nn*(1-p-q)^2 + .175*(nc+cn)*q*(1-p-q) + .325*(ns+sn)*p*(1-p-q) + .40*cc*q^2 + .545*(sc+cs)*p*q + .42*ss*p^2)/norm(p, q)
|
|
def F(Z):
(p, q) = Z
return (Pt(p, q), Qt(p, q))
|
|
p0 = 0.011
q0 = 0.027
t = range(101)
Z = (p0, q0)
pt = [p0]
qt = [q0]
for T in t:
Z = F(Z)
pt.append(Z[0])
qt.append(Z[1])
|
|
line(zip(t, pt)) + line(zip(t, qt), color = 'red')
|
|
equilibrium = (pt[-1], qt[-1])
equilibrium
(0.0208332371765571, 0.0513933720444454) (0.0208332371765571, 0.0513933720444454) |
J = jacobian([Pt(p, q), Qt(p, q)], [p, q])(p = equilibrium[0], q = equilibrium[1])
J
[0.203318214372026 0.112225078298775] [0.490862107783925 0.278908686176197] [0.203318214372026 0.112225078298775] [0.490862107783925 0.278908686176197] |
map(n, J.eigenvalues())
[0.00338352319575858, 0.478843376838240] [0.00338352319575858, 0.478843376838240] |
pt
[0.0110000000000000, 0.0161355834145870, 0.0185863573678157, 0.0197579273190266, 0.0203184690126116, 0.0205867753055159, 0.0207152277628683, 0.0207767308070029, 0.0208061798556721, 0.0208202810447708, 0.0208270332386326, 0.0208302664665485, 0.0208318146727933, 0.0208325560202911, 0.0208329110094453, 0.0208330809936083, 0.0208331623893893, 0.0208332013652178, 0.0208332200285346, 0.0208332289653402, 0.0208332332446703, 0.0208332352937992, 0.0208332362750110, 0.0208332367448577, 0.0208332369698407, 0.0208332370775724, 0.0208332371291589, 0.0208332371538608, 0.0208332371656892, 0.0208332371713531, 0.0208332371740652, 0.0208332371753639, 0.0208332371759858, 0.0208332371762836, 0.0208332371764261, 0.0208332371764944, 0.0208332371765271, 0.0208332371765428, 0.0208332371765503, 0.0208332371765538, 0.0208332371765556, 0.0208332371765564, 0.0208332371765568, 0.0208332371765570, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571] [0.0110000000000000, 0.0161355834145870, 0.0185863573678157, 0.0197579273190266, 0.0203184690126116, 0.0205867753055159, 0.0207152277628683, 0.0207767308070029, 0.0208061798556721, 0.0208202810447708, 0.0208270332386326, 0.0208302664665485, 0.0208318146727933, 0.0208325560202911, 0.0208329110094453, 0.0208330809936083, 0.0208331623893893, 0.0208332013652178, 0.0208332200285346, 0.0208332289653402, 0.0208332332446703, 0.0208332352937992, 0.0208332362750110, 0.0208332367448577, 0.0208332369698407, 0.0208332370775724, 0.0208332371291589, 0.0208332371538608, 0.0208332371656892, 0.0208332371713531, 0.0208332371740652, 0.0208332371753639, 0.0208332371759858, 0.0208332371762836, 0.0208332371764261, 0.0208332371764944, 0.0208332371765271, 0.0208332371765428, 0.0208332371765503, 0.0208332371765538, 0.0208332371765556, 0.0208332371765564, 0.0208332371765568, 0.0208332371765570, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571, 0.0208332371765571] |
qt
[0.0270000000000000, 0.0398027440308597, 0.0458639215325515, 0.0487503649405818, 0.0501288709345327, 0.0507881229132987, 0.0511036095646767, 0.0512546342828259, 0.0513269413858504, 0.0513615628513294, 0.0513781405806865, 0.0513860785950553, 0.0513898796327968, 0.0513916997281569, 0.0513925712673013, 0.0513929885977130, 0.0513931884335399, 0.0513932841235845, 0.0513933299441246, 0.0513933518849859, 0.0513933623912218, 0.0513933674220632, 0.0513933698310483, 0.0513933709845749, 0.0513933715369334, 0.0513933718014267, 0.0513933719280775, 0.0513933719887234, 0.0513933720177633, 0.0513933720316689, 0.0513933720383274, 0.0513933720415158, 0.0513933720430426, 0.0513933720437737, 0.0513933720441238, 0.0513933720442914, 0.0513933720443717, 0.0513933720444101, 0.0513933720444285, 0.0513933720444373, 0.0513933720444415, 0.0513933720444436, 0.0513933720444445, 0.0513933720444450, 0.0513933720444452, 0.0513933720444453, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454] [0.0270000000000000, 0.0398027440308597, 0.0458639215325515, 0.0487503649405818, 0.0501288709345327, 0.0507881229132987, 0.0511036095646767, 0.0512546342828259, 0.0513269413858504, 0.0513615628513294, 0.0513781405806865, 0.0513860785950553, 0.0513898796327968, 0.0513916997281569, 0.0513925712673013, 0.0513929885977130, 0.0513931884335399, 0.0513932841235845, 0.0513933299441246, 0.0513933518849859, 0.0513933623912218, 0.0513933674220632, 0.0513933698310483, 0.0513933709845749, 0.0513933715369334, 0.0513933718014267, 0.0513933719280775, 0.0513933719887234, 0.0513933720177633, 0.0513933720316689, 0.0513933720383274, 0.0513933720415158, 0.0513933720430426, 0.0513933720437737, 0.0513933720441238, 0.0513933720442914, 0.0513933720443717, 0.0513933720444101, 0.0513933720444285, 0.0513933720444373, 0.0513933720444415, 0.0513933720444436, 0.0513933720444445, 0.0513933720444450, 0.0513933720444452, 0.0513933720444453, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454, 0.0513933720444454] |
derivative(Pt, q)
(p, q) |--> -(-0.212000000000000*p - 0.212000000000000*q - 0.148000000000000)*(0.0102740000000000*(p + q - 1)^2 - 0.200200000000000*(p + q - 1)*p - 0.120400000000000*(p + q - 1)*q + 0.190000000000000*p^2 + 0.300900000000000*p*q + 0.108800000000000*q^2)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2)^2 + (0.000848000000000015*p - 0.00265199999999999*q + 0.0998520000000000)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2) (p, q) |--> -(-0.212000000000000*p - 0.212000000000000*q - 0.148000000000000)*(0.0102740000000000*(p + q - 1)^2 - 0.200200000000000*(p + q - 1)*p - 0.120400000000000*(p + q - 1)*q + 0.190000000000000*p^2 + 0.300900000000000*p*q + 0.108800000000000*q^2)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2)^2 + (0.000848000000000015*p - 0.00265199999999999*q + 0.0998520000000000)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2) |
derivative(Qt, q)
(p, q) |--> -(-0.212000000000000*p - 0.212000000000000*q - 0.148000000000000)*(0.0252180000000000*(p + q - 1)^2 - 0.500500000000000*(p + q - 1)*p - 0.301000000000000*(p + q - 1)*q + 0.210000000000000*p^2 + 0.643100000000000*p*q + 0.272000000000000*q^2)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2)^2 + (-0.107964000000000*p - 0.00756399999999990*q + 0.250564000000000)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2) (p, q) |--> -(-0.212000000000000*p - 0.212000000000000*q - 0.148000000000000)*(0.0252180000000000*(p + q - 1)^2 - 0.500500000000000*(p + q - 1)*p - 0.301000000000000*(p + q - 1)*q + 0.210000000000000*p^2 + 0.643100000000000*p*q + 0.272000000000000*q^2)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2)^2 + (-0.107964000000000*p - 0.00756399999999990*q + 0.250564000000000)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2) |
derivative(Qt, p)
(p, q) |--> (-0.530564000000000*p - 0.107964000000000*q + 0.450064000000000)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2) - (-0.212000000000000*p - 0.212000000000000*q - 0.328000000000000)*(0.0252180000000000*(p + q - 1)^2 - 0.500500000000000*(p + q - 1)*p - 0.301000000000000*(p + q - 1)*q + 0.210000000000000*p^2 + 0.643100000000000*p*q + 0.272000000000000*q^2)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2)^2 (p, q) |--> (-0.530564000000000*p - 0.107964000000000*q + 0.450064000000000)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2) - (-0.212000000000000*p - 0.212000000000000*q - 0.328000000000000)*(0.0252180000000000*(p + q - 1)^2 - 0.500500000000000*(p + q - 1)*p - 0.301000000000000*(p + q - 1)*q + 0.210000000000000*p^2 + 0.643100000000000*p*q + 0.272000000000000*q^2)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2)^2 |
derivative(Pt, p)
(p, q) |--> -(-0.212000000000000*p - 0.212000000000000*q - 0.328000000000000)*(0.0102740000000000*(p + q - 1)^2 - 0.200200000000000*(p + q - 1)*p - 0.120400000000000*(p + q - 1)*q + 0.190000000000000*p^2 + 0.300900000000000*p*q + 0.108800000000000*q^2)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2)^2 + (0.000147999999999981*p + 0.000848000000000015*q + 0.179652000000000)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2) (p, q) |--> -(-0.212000000000000*p - 0.212000000000000*q - 0.328000000000000)*(0.0102740000000000*(p + q - 1)^2 - 0.200200000000000*(p + q - 1)*p - 0.120400000000000*(p + q - 1)*q + 0.190000000000000*p^2 + 0.300900000000000*p*q + 0.108800000000000*q^2)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2)^2 + (0.000147999999999981*p + 0.000848000000000015*q + 0.179652000000000)/(0.934000000000000*(p + q - 1)^2 - 1.54000000000000*(p + q - 1)*p - 1.72000000000000*(p + q - 1)*q + 0.500000000000000*p^2 + 1.18000000000000*p*q + 0.680000000000000*q^2) |
var('p` q` p`` q``')
(p`, q`, p``, q``) (p`, q`, p``, q``) |