# How to integrate over a set of data?

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
$\int_{-\infty}^{\infty}dx_0\int_{-\infty}^{\infty}dp_0e^{-{\beta}H(x_0,p_0)}x_0x(t)$.
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?

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:

TMAX=10;
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, -∞,∞}]
]


Edit:

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, -∞,∞}]

• Thank you. This works great. I didn't know I could get a numerical solution with parameters before. Apr 26, 2017 at 4:36