# Simulating large nontrivial SDE systems: algorithms, reliability and speed

**URL:** <https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082>\
**Category:** Modelling & Simulations\
**Tags:** question, diffeq, sde\
**Created:** [June 30, 2018, 2:40pm UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082 "2018-06-30T14:40:37Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![lappet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lappet/32/4566_2.png) [@lappet](https://discourse.julialang.org/u/lappet)\
**Post date:** [June 30, 2018, 2:40pm UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082/1 "2018-06-30T14:40:37Z")

</div>

In short, I am interested in the simulation of large systems of nontrivial **stochastic differential equations** as they arise when applying phase-space methods to **interacting quantum systems**.

As an intro example: when one uses the so-called “positive P representation” on the two-site Bose-Hubbard model, this would result in **eight nonlinear equations** for the deterministic part, along with **four independent Wiener processes** that then enter the system **as multiplicative noises** , i.e `noise_rate_prototype=zeros(8, 4)` .

Note: I am _not_ asking how to solve this specific system.

Rather, I would like to have some input on the more general questions:

1. What algorithms are suggested for systems that have both, **nonlinearity** _and_ **multiplicative noise?** (For the two-site Bose-Hubbard, I have tested `ImplicitEM()`, with satisfactory results in ‘unproblematic’ parameter regimes.)

2. Numerically, how large can I make such systems and still get **a reliable solution in finite time?** For the combination of Bose-Hubbard and positive P, for instance, the degrees of freedom will have a factor of four w.r.t. the number of sites, and I will have half as many independent Wiener processes. Say I have 128 sites, that will make **512 degrees of freedom and 256 Wieners** , combined in a coupled system of equations of the type described above. My naive intuition suggests that, among other things, one will have to create ever larger ensembles when the number of noises goes up. The nonlinearity will add further complications, obviously. Hence the issue of **speed vs. accuracy** seems to be of great concern here, too.

Thanks in advance for any help!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 1, 2018, 3:27am UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082/2 "2018-07-01T03:27:22Z")

</div>

Great question with lots of answers to explore.

## General Overview

The current state of where we are at is well-summarized in this preprint:

> **[Stability-Optimized High Order Methods and Stiffness Detection for Pathwise...](https://arxiv.org/abs/1804.04344)**
>
> Stochastic differential equations (SDE) often exhibit large random transitions. This property, which we denote as pathwise stiffness, causes transient bursts of stiffness which limit the allowed step size for common fixed time step explicit and...

Not all of the algorithms are in there, because for example you’ll notice that there are some algorithms here which have no source:

[http://docs.juliadiffeq.org/latest/solvers/sde\_solve.html#Stiff-Methods-1](http://docs.juliadiffeq.org/latest/solvers/sde_solve.html#Stiff-Methods-1)

I will get those written up, but they are adaptive and for the large noise case since methods like `ImplicitEM()` are still largely unstable when the noise term is sufficient large (since implictness on SDE solvers is semi-implicit only in the drift term due to the non-invertability of standard normal random variables, though there are some truncation schemes that can be done and this is quite a rabbit hole I would like to go down but not right now).

But okay, so the state is basically that if you have an additive noise SDE there’s a ton of stuff we can do, and that paper shows how to convert an affine noise (or multiplicative noise) SDE into an additive noise SDE (under certain conditions for invertability of the noise term). If it’s a stiff diagonal noise SDE, in many cases `SOSRI` is doing better than the (semi) implicit algorithms which is showing that there may be a deficiency in their design, something that has been noted before:

[https://www.newton.ac.uk/files/preprints/ni07050.pdf](https://www.newton.ac.uk/files/preprints/ni07050.pdf)

We have plans (and students!) to try out alternatives like the SROCK methods, track it here:

[https://github.com/JuliaDiffEq/StochasticDiffEq.jl/issues/73](https://github.com/JuliaDiffEq/StochasticDiffEq.jl/issues/73)

## What it means for you

Well there’s a few more things to mention. When you say four independent Wiener processes, do you mean real-valued Wiener processes or complex-valued Wiener processes? If they are complex valued, then your SDE isn’t diagonal but actually commutative (since multiplication mixes the real and complex valued terms) and so the diagonal noise methods will have order loss to 1.0. This is another reason to say that by far the best technique will be to convert to an additive system using the Lamperti transform mentioned in the paper, and then applying `SOSRA` or `SKenCarp`.

One thing we do need to add to our implementation though is the ability to tell it how many random numbers you want. I haven’t seen this factor into profiles though, and the diagonal/additive noise methods are so much more efficient than the general SDE methods like `ImplicitEM` that I would be heavily surprised if using `noise_rate_prototype=zeros(8, 4)` is actually helpful here. On the topic of profiling…

## Please donate an example problem

Most of the testing and development has centered around a set of 20 highly stiff biological SDEs because someone gave this nice problem to me:

[https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/blob/master/src/sde\_premade\_problems.jl#L172-L318](https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/blob/master/src/sde_premade_problems.jl#L172-L318)

Please feel free to donate an implementation of this problem for us to benchmark and profile on. That would be very helpful for algorithm development since for SDEs it’s not really known what all of the properties of equations that are important for real-world performance, so testing directly on problems related to the application can be a good idea.

Also, if you cannot convert it to an additive noise problem due to invertability issues and it’s a complex Wiener process, that would be a great real-world commutative noise example and I would really like to have an implementation of that for testing some algorithms in development!

---

<div class="post-metadata">

**Author:** ![lappet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lappet/32/4566_2.png) [@lappet](https://discourse.julialang.org/u/lappet)\
**Post date:** [July 2, 2018, 8:21am UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082/3 "2018-07-02T08:21:54Z")

</div>

Thank you for the detailed answer! I will post the problem here after I have digested your input!

---

<div class="post-metadata">

**Author:** ![lappet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lappet/32/4566_2.png) [@lappet](https://discourse.julialang.org/u/lappet)\
**Post date:** [July 2, 2018, 5:15pm UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082/4 "2018-07-02T17:15:24Z")

</div>

Here is the example implementation of the problem. When I said “four independent Wiener processes”, I actually did mean real increments, although I realize now that I probably should not have called the noise “multiplicative”… It is, in some sense, as you can see below, but seemingly not in the convenient “diagonal” form that would lend itself immediately to a Lamperti transform, as far as I can see. This was probably expected since simulating interacting quantum systems should not be that easy.

Concerning the parameters I have set below, it should be noted that they are very small: one certainly wants the solution of these equations for either large U or N, and possibly also for cases where both are not small. Furthermore, the deterministic system has an instability in terms of the dimensionless parameter \Lambda = NU/2J, as well as a so-called “self-trapped” regime for large initial values of z. In this sense, the simple two-site model is already very rich and a good playground to start off from.

As I said, I already have obtained some results for the parameters below with `ImplicitEM()`, but I certainly did have problems stretching that approach to the more interesting parameter regimes…

```julia
# some arbitrary evolution time
t_max = 1.

# hopping
J = -2*pi*1.0/2.

# interaction (this is a very weak value, much larger ones definitely desirable)
U = 2*pi*0.005

# (very small) number of particles (should be approx. conserved)
N = 10.0

# deterministic part (linear terms prop. to J, nonlinear terms prop. to U)
function sys(du,u,p,t)
  du[1] = J*u[6] + U*u[1]^2*u[7] + 2*U*u[1]*u[3]*u[5] - U*u[5]^2*u[7]
  du[2] = J*u[5] + U*u[2]^2*u[8] + 2*U*u[2]*u[4]*u[6] - U*u[6]^2*u[8]
  du[3] = -J*u[8] - 2*U*u[1]*u[3]*u[7] - U*u[3]^2*u[5] + U*u[5]*u[7]^2
  du[4] = -J*u[7] - 2*U*u[2]*u[4]*u[8] - U*u[4]^2*u[6] + U*u[6]*u[8]^2
  du[5] = -J*u[2] - U*u[1]^2*u[3] + 2*U*u[1]*u[5]*u[7] + U*u[3]*u[5]^2
  du[6] = -J*u[1] - U*u[2]^2*u[4] + 2*U*u[2]*u[6]*u[8] + U*u[4]*u[6]^2
  du[7] = J*u[4] + U*u[1]*u[3]^2 - U*u[1]*u[7]^2 - 2*U*u[3]*u[5]*u[7]
  du[8] = J*u[3] + U*u[2]*u[4]^2 - U*u[2]*u[8]^2 - 2*U*u[4]*u[6]*u[8]
end

# noise strength (related to interaction)
sigma = sqrt(U)/2

# nondiagonal (?) multiplicative noise: the Wiener increments
# for du[1, 1] and du[5, 1] are the same, and similarly for
# the remaining three pairs (du[2, 2] and du[6, 2] etc.)
function σ_sys(du,u,p,t)
    du[1, 1] = -sigma*(u[1] + u[5])
    du[2, 2] = -sigma*(u[2] + u[6])
    du[3, 3] = sigma*(u[3] - u[7])
    du[4, 4] = sigma*(u[4] - u[8])
    du[5, 1] = sigma*(u[1] - u[5])
    du[6, 2] = sigma*(u[2] - u[6])
    du[7, 3] = sigma*(u[3] + u[7])
    du[8, 4] = sigma*(u[4] + u[8])
end

# so-called population imbalance z = (N_1 - N_2)/(N_1 + N_2)
z = 0.1

# initial population derived from z
N_1 = 0.5*N*(1+z)
N_2 = 0.5*N*(1-z)

# SDE problem where initial phases are set to zero, hence real initial values
prob_sde = SDEProblem(sys, σ_sys, [sqrt(N_1), sqrt(N_2), sqrt(N_1), sqrt(N_2), 0., 0., 0., 0.], (0.0, t_max), noise_rate_prototype=zeros(8, 4))

```

---

<div class="post-metadata">

**Author:** ![lappet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lappet/32/4566_2.png) [@lappet](https://discourse.julialang.org/u/lappet)\
**Post date:** [July 2, 2018, 5:25pm UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082/5 "2018-07-02T17:25:29Z")

</div>

Let me remark a couple of other things.

Concerning the people who invented the mapping from quantum to classical description I am using above, they have provided an algorithm [here](https://www.sciencedirect.com/science/article/pii/S0021999196956387). There is also a software package related to that called `xmds`. (I can only post two links for a start?) The algorithm is still in use in [more recent work](https://arxiv.org/pdf/1412.5579.pdf).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 2, 2018, 6:25pm UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082/6 "2018-07-02T18:25:30Z")

</div>

[Robust Algorithms for Solving Stochastic Partial Differential Equations - ScienceDirect](https://www.sciencedirect.com/science/article/pii/S0021999196956387) is fine. The algorithm they build in to 2.9 is an implicit form of `RKMilCommute(interpretation=:Stratanovich)`. Does your model satisfy the commutivity relation? It wouldn’t be hard to implement the implicit form (it’s mostly pasting together `ImplicitRKMil` and `RKMilCommute` and adding tests) if you need it.

[http://docs.juliadiffeq.org/latest/solvers/sde\_solve.html#Special-Noise-Forms-1](http://docs.juliadiffeq.org/latest/solvers/sde_solve.html#Special-Noise-Forms-1)

If it does then that’s likely the method you want to be using. Or for large `sigma` you might want to use `ISSEM` because the semi-implicit methods (like `ImplicitEM` and the `ImplicitRKMilCommute` from the paper) are not stable for large noise eigenvalues, while step-splitting corrects for this. I think what you really want is an `ISSRKMilCommute`.

So please do check and if it does satisfy the commutivity equation I can add these algorithms since this is really what you need and I want to start researching better methods for this case anyways.

The rest of the paper is just the spatial discretization which is part of the SDE definition so I’ll leave that for you.

---

<div class="post-metadata">

**Author:** ![lappet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lappet/32/4566_2.png) [@lappet](https://discourse.julialang.org/u/lappet)\
**Post date:** [July 3, 2018, 7:06am UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082/7 "2018-07-03T07:06:16Z")

</div>

Okay, so it was quite easy to see that the noise _is_ indeed commutative. I get that it makes a lot of sense then to exploit this property algorithmically.

Thanks again for your great support!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 3, 2018, 7:30am UTC](https://discourse.julialang.org/t/simulating-large-nontrivial-sde-systems-algorithms-reliability-and-speed/12082/8 "2018-07-03T07:30:43Z")

</div>

> [@lappet](#):
>
> Okay, so it was quite easy to see that the noise _is_ indeed commutative. I get that it makes a lot of sense then to exploit this property algorithmically.

Yup, it means the Wiktersson approximation of the iterated integrals decomposes and then is much simpler to compute.