Solving Stochastic Differential Equations with Same Random Numbers


One Line Question: How do I use Julia to solve a system of Stochastic Differential Equations were there are correlated terms of noise (or “the same” noise appears in several places)

Biological Background: Lets say I have two components, X and Y. Now we have a reaction X <-> Y. I want to model this using a standard approach. I also want to include the stochastics of the system using the standard Langevin approach. I get (Stochastic) Differential Equations:
d[X]/dt = r[Y] - q[X] + sqrt(r[Y])*N1(0,1) - sqrt(q[X])*N2(0,1)
d[Y]/dt = q[Y] - r[X] + sqrt(q[X])*N2(0,1) - sqrt(r[Y])*N1(0,1)
Here N1 and N2 are normally distributed Noise with mean 0 and std 1. What is important is that the same N1 and N2 appears in both equations (which in this case has the important implication that the sum [X] + [Y] is preserved.

I (think) I understand how to use non-diagonal noise and the stochastic differential equation solver to solve this problem, but then all the 4 instances of noise would all be uncorrelated
( ). However, I would I do to solve the problem were the same noise appears in several places?

(The example above is not the one I am exactly interested in, but a simple example of the Chemical Langevin Equations were I do not know how to solve it in Julia. However, if I understand this example I should be able to do int in a more complicated case as well).

Thanks in advance


Note that I took the liberty of changing the title because correlated noise is different so hopefully it helps people find this in searches better in the future

This is what non-diagonal noise is for. Let’s look at the example from the documentation:

f(t,u,du) = du .= 1.01u
function g(t,u,du)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))

This is the equation

du = f(t,u)dt + g(t,u)*dW

where W is a vector of 4 noise terms. This means if you look at the matrix multiplication, g(t,u)*dW, then that means that the SDEs are:

du1 = f1(t,u)dt + g11(t,u)*dW1 + g12(t,u)*dW2 + g13(t,u)*dW3 + g14(t,u)*dW4
du2 = f1(t,u)dt + g21(t,u)*dW1 + g22(t,u)*dW2 + g23(t,u)*dW3 + g24(t,u)*dW4

As you can see by writing it out, g[1,1] and g[2,1] are two terms which have the same random number applied to them. Another way to explain this is if I simplify this example a bit:

f(t,u,du) = du .= 1.01u
function g(t,u,du)
  du[1,1] = 0.3u[1]
  du[2,1] = 1.2u[1]
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,1))

then by the same logic this is an SDE with a single Brownian motion and both SDEs get the same random number.

This blog post shows exactly your same problem (chemical Langivan equations and non-diagonal noise to preserve sums). Hopefully this helps as well.

Hopefully together these explain it. But since it the docs didn’t work the first time, please let us know how to improve the documentation so that this is more clear!


Thanks for your reply, I understand.

Now when I read the documentation I understand it well, but I am unsure why I was not able to at the time. Partially the reason should have been that I was unfamiliar with Diagonal/Non-Diagonal Noise. Maybe when I looked at the two sections the obvious difference I saw was that in one you had several noise terms in each equation and in one not, and missed that in the Non-Diagonal case the various noise terms were the same.


@Gaussia since this explanation worked I updated the docs to say something similar. Thanks for reporting this since hopefully this clears it up for more people!