Issue: Correctly Implementing the Spherical Constraint
I am working on simulating the spherical p=2 spin model using the DifferentialEquations.jl package in Julia. The model consists of N coupled stochastic differential equations (SDEs):
\dot{x}_i(t) = -\lambda(t)x_i(t) + \sum_{j=1}^{N} J_{ij}x_j(t)+\sqrt{2 D}\eta_i(t)
where:
- x_i(t) is the state of the i-th spin,
- J_{ij} is the sparse coupling matrix,
- \eta_i(t) is a Gaussian white noise (zero mean, variance 1, independent across sites and times),
- \lambda(t) is a Lagrange multiplier that enforces the spherical constraint, given by:
\lambda(t) = \frac{1}{N} \sum_{i=1}^N x_i(t) \left( \eta_i(t) + \sum_{j \in \partial i} J_{ij} x_j(t) \right)
Current Approach
I have implemented a simulation without enforcing the spherical constraint and am now trying to correctly update \lambda(t) at each time step. Initially, I considered using a callback, but I now prefer directly modifying the integrator’s state inside the solver.
One challenge I am facing is that the noise realization dW is generated internally by the solve function, so I cannot directly pass it to a function f! that updates the deterministic part of the equation in-place. What is the best way to enforce the spherical constraint while ensuring that the correct noise realization is used at each step? How can I compute \lambda(t) without interfering with the solver’s internal noise generation?