I would like to define an Hamiltonian dynamics that I will integrate later using scipy.integrate. The Hamiltonian is computed with a symbolic expression and then I need to make substitution from the symbolic variables to the numerical variables.
N=2
var('q1 q2')
var('p1 p2')
qs = [q1,q2]
ps = [p1,p2]
def dynq(t,z):
H = p1*(q1^2+1/3*q2^2) + p2*(cos(q1)-2*sin(q2*q1)^2)
jacHp = jacobian(H,tuple(ps))
dqdt = list(jacHp)
dqdt = list(dqdt[0])
dqdt[0].subs({qs[i]:z[i] for i in range(0,N)})
dqdt[1].subs({qs[i]:z[i] for i in range(0,N)})
print(dqdt)
return dqdt
I call dynq(0.,[1.,0.5])
and the output is
[q1^2, -2*sin(q1*q2)^2 + cos(q1)]
so q1,q2,p1,p2 are still in the expression and are not replaced by the numerical values of the list $z$.