Numerical solution of Parisi equations for the SK spin-glass model

Question

I am trying to reproduce the numerical solution of the so-called Parisi equations for the Sherrington-Kirkpatrick (SK) model. There are at least two main methods to solve these equations, and unfortunately I am confused by both.

First, the equations themselves. They are:

$$ q(x) = \int_{-\infty}^{\infty} dy \, P(x,y) m(x,y)^2 \,,$$

$$ \dot{m}(x,y) = - \frac{\dot{q}(x)}{2} \left[ m''(x,y) + 2 \beta x m(x,y) m'(x,y) \right] \,,$$

$$ \dot{P}(x,y) = \frac{\dot{q}(x)}{2} \left[ P''(x,y) - 2 \beta x \left[ m(x,y) P(x,y) \right]' \right] \,,$$

and the initial conditions are:

$$ m(1,y) = \tanh(\beta y), \qquad P(0,y) = \delta(y) \,.$$

Primes denote $y$ derivatives and a dot denotes an $x$ derivative. The domain is $x \in [0,1]$, $y \in [-\infty, \infty]$. In any numerical attempt, the range of $y$ is approximated as finite, $y \in [-y_{\text{max}}, y_{\text{max}}]$. 

In each of the ways to solve these equations that I know about, one proceeds by first supplying a guess for $q(x)$, then solving the $m$ equation, then the $P$ equation. Then $q$ is updated according to the first equation. These steps are iterated until convergence, which typically takes hundreds of cycles .

The two main methods I know of are an integral kernel method by Nemoto and a pseudo-spectral method by Crisanti and Rizzo. In the Nemoto method kernels are used, and these require an initial guess for not only $q(x)$ but also $m(x,y), P(x,y)$. I have no idea what a good guess would be, and I have tried many and failed to get convergence. Presumably the guess should agree with the initial conditions for $m, P$ for $x=1, 0$ respectively. The Crisanti and Rizzo method takes a Fourier transform in the $y$-direction, and then numerically integrates in the $x$-direction using an Adams-Bashforth scheme. Implicit in their approach is that they are setting $P, m$ to have Dirichlet boundary conditions in the $y$-direction, although this is inconsistent with the initial condition for $m$. I have tried to reproduce their method using using both Mathematica's NDSolve, and by adapting this pseudo-spectral code by Trefethen for the KdV equation, but I cannot get either code to converge. From the shape of the curves it is clear I am not handling the boundary conditions properly.

I would greatly appreciate any guidance, as well as any references that might shed light or offer a more easy-to-understand numerical method to solve these equations.

Answer 1

I have a partial answer and a question. 

I am also trying to compute a numerical solution to the Parisi formula for the SK model (with the hope of simulating Montanari's algorithm; is this your motivation as well?). While I can roughly solve for any fixed $\beta$ (see below if you're curious for how to do that), I do not know how to take a limit as $\beta$ goes to infinity as one must do to find the ground states. This seems numerically difficult. I was wondering if you had any thoughts on this / have succeeded in the simulation since posting the question?

Anyway, I think a simple way to solve for the Parisi measure at a fixed temperature is to break the problem into two parts: for any given measure, one can use the Matlab "pdepe" solver (applicable since the PDE is hyperbolic and just in 1 space variable). This allows us to evaluate the Parisi functional for any measure, so then we can use stochastic optimization algorithms, e.g. simulated annealing / genetic algorithms, to search over the space of feasible measures.