"""
War - Lanchester's model

code for slide talk by David Joyner, 2011-04-22
license: public domain
"""

A = matrix([[0, -1], [-25, 0]]) 
A.eigenspaces_right()
"""
[ (5, Vector space of degree 2 and dimension 1 over 
Rational Field User basis matrix: [ 1 -5]), 
  (-5, Vector space of degree 2 and dimension 1 over 
Rational Field User basis matrix: [1 5]) ]
"""

A.eigenvectors_right()

c1,c2 = var("c1,c2") 
solve([c1+c2==27, -5*c1+5*c2==33],[c1,c2]) 
# [[c1 == (51/5), c2 == (84/5)]]


"""

t = var('t')    
x = function('x', t)
y = function("y", t)
z = function("z", t)
de1 = diff(x,t) == -x - 2*z
de2 = diff(y,t) == -2*x - z
de3 = diff(z,t) == -x - 2*y
desolve_system([de1, de2, de3], [x,y,z])

desolve_system([de1, de2, de3], [x,y,z], ics=[0,1,1,1])

# x vs y vs z- tie

desolve_system([de1, de2, de3], [x,y,z], ics=[0,1,0,1])

Px = plot(-4/17*(4*sin(t) - cos(t))*e^t + 13/17*e^(-3*t),(t,0,1),rgbcolor=(1,0,0),legend_label='x')
Py = plot(1/17*(sin(t) - 13*cos(t))*e^t + 13/17*e^(-3*t),(t,0,1),rgbcolor=(0,1,0),legend_label='y')
Pz = plot(2/17*(9*sin(t) + 2*cos(t))*e^t + 13/17*e^(-3*t),(t,0,1),rgbcolor=(0,0,1),legend_label='z')
show(Px+Py+Pz)
#z vs x - z wins

desolve_system([de1, de2, de3], [x,y,z], ics=[0,1,1,0])

Px = plot(2/17*(5*sin(t) + 3*cos(t))*e^t + 11/17*e^(-3*t),(t,0,1),rgbcolor=(1,0,0),legend_label='x')
Py = plot(-1/17*(7*sin(t) - 6*cos(t))*e^t + 11/17*e^(-3*t),(t,0,1),rgbcolor=(0,1,0),legend_label='y')
Pz = plot(-1/17*(7*sin(t) + 11*cos(t))*e^t + 11/17*e^(-3*t),(t,0,1),rgbcolor=(0,0,1),legend_label='z')
show(Px+Py+Pz)
#y vs x - x wins

desolve_system([de1, de2, de3], [x,y,z], ics=[0,0,1,1])
Px = plot(2/17*(3*sin(t) - 5*cos(t))*e^t + 10/17*e^(-3*t),(t,0,1),rgbcolor=(1,0,0),legend_label='x')
Py = plot(1/17*(6*sin(t) + 7*cos(t))*e^t + 10/17*e^(-3*t),(t,0,1),rgbcolor=(0,1,0),legend_label='y')
Pz = plot(-1/17*(11*sin(t) - 7*cos(t))*e^t + 10/17*e^(-3*t),(t,0,1),rgbcolor=(0,0,1),legend_label='z')
show(Px+Py+Pz)
#y vs z - y wins

"""


"""
Love - Romeo and Juliet model
"""

s,t = var("s,t")
A = matrix(SR, [[s,-5,4],[1,s-2,6]])
B = A.echelon_form() 
Rs = B[0][2]; Js = B[1][2]
rt = Rs.inverse_laplace(s,t); rt 
# (13*sin(2*t) + 4*cos(2*t))*e^t 
jt = Js.inverse_laplace(s,t); jt 
# (sin(2*t) + 6*cos(2*t))*e^t 
parametric_plot((rt,jt),(t,0,1.5))


t = var("t") 
r = function("r",t) 
j = function("j",t) 
de1 = diff(r,t) == 5*j 
de2 = diff(j,t) == -r+(1/5)*j 
soln = desolve_system([de1, de2], [r,j], ics=[0,4,6]) 
rt = soln[0].rhs(); jt = soln[1].rhs() 
parametric_plot((rt,jt),(t,0,5))



"""
Zombies 

"""

from sage.calculus.desolvers import desolve_system_rk4 sage: t,s,z,r = var(’t,s,z,r’) 
a,b,zeta,d,B = 0.005,0.0095,0.0001,0.0001,0.0 
P = desolve_system_rk4([B-b*s*z-d*s,b*s*z-zeta*r-a*s*z,d*s+a*s*z-zeta*r],[s,z,r],ics=[0,11,10,5] ,ivar=t,end_points=30) 
Ps = list_plot([[t,s] for t,s,z,r in P],plotjoined=True,legend_label=’People’) 
Pz = list_plot([[t,z] for t,s,z,r in P],plotjoined=True,rgbcolor=’red’,legend_label=’Zombies’)
Pr = list_plot([[t,r] for t,s,z,r in P],plotjoined=True,rgbcolor=’black’,legend_label=’Deceased’)
show(Ps+Pz+Pr)