Solving Stochastic Differential Equations with Same Random Numbers

Specific Domains Modelling & Simulations
,
1

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/dt = r[Y] - q + sqrt(r[Y])*N1(0,1) - sqrt(q)*N2(0,1)
d[Y]/dt = q[Y] - r + sqrt(q)*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 + [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
( http://docs.juliadiffeq.org/latest/tutorials/sde_example.html ). 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

2

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:

http://docs.juliadiffeq.org/latest/tutorials/sde_example.html#Example-4:-Systems-of-SDEs-with-Non-Diagonal-Noise-1

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]
end
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]
end
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.

http://mbesancon.github.io/post/2017-12-14-diffeq-julia/

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!

3

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.

4

@Torkel 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!