The Hamiltonian of a system is $H = 0.5p^2 + 0.5x^2 + 0.1x^4$. I have obtained the equations of motion. To solve them (numerically) I need to specify $x_0$ and $p_0$. But then the question asks me to compute
Since $x(t)$ depends on the values of $x_0$ and $p_0$ that I specify, I have to integrate then over a list of data. How do I do this?


1 Answer 1


I'm not quite sure that this is what you're trying to do, but if I do then here is my suggestion:

I assume the equations of motion are $\dot x=p$ and $\dot p=-x-\frac{2}{5}x^3$. With these equations, the value of $e^{-\beta H}x(0)x(T)$ can be obtained using ParametricNDSolveValue as follows:

fun = ParametricNDSolveValue[{x'[t] == p[t], 
   p'[t] == -x[t] - 2/5 x[t]^3, x[0] == x0, p[0] == p0}, 
   Exp[-β (1/2 (p[T]^2 + x[T]^2) + 1/10 x[T]^4)] x[T], 
   {t, 0, TMAX}, {x0, p0}];

To get the value of $e^{-\beta H}x(0)x(T)$ for $x_0=1,p_0=2$ and $t=3$ you simply evaluate fun[1,2,3]. Note, however, that this is done by numerically solving the equations from t=0 to TMAX so fun can not be evaluated for t>TMAX (technically it does return a value, but it's an extrapolation and you shouldn't use that).

Then you can simply integrate. For example, to get the answer at $t=2$ and $\beta=1$ you can run

Block[{β = 1}, 
   NIntegrate[fun[x0, p0, 2], {x0, -∞,∞}, {p0, -∞,∞}]


Now I see that you want the Hamiltonian to be evaluated at $x_0,p_0$, so you can use a simpler code. You replace fun by

fun = ParametricNDSolveValue[{x'[t] == p[t], 
  p'[t] == -x[t] - 2/5 x[t]^3, x[0] == x0, p[0] == p0}, 
  x[T], {t, 0, TMAX}, {x0, p0, T}]

and you can set up a function to calculate the correlation like this:

corr[β_, t_]:= NIntegrate[
    Exp[-β (1/2 (p0^2 + x0^2) + 1/10 x0^4)] x0 fun[x0, p0,t], 
    {x0, -∞,∞}, {p0, -∞,∞}]
  • $\begingroup$ Thank you. This works great. I didn't know I could get a numerical solution with parameters before. $\endgroup$
    – user21417
    Apr 26, 2017 at 4:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.