Some inconsistencies in an ensemble of stochastic differential equations

General Usage
,
1

I am solving a stochastic differential equation using EnsembleProblem.
When I then simply plot a summary of the solution, I get a reasonable result.
However, when I calculate the mean myself using timeseries_point_meanvar I find that the same solution is unstable.

Here is the ensemble summary
image

And here is the calculated mean with timeseries_point_meanvar
image

The code I am running for the 1st plot is:

using DifferentialEquations
using Plots

function rayleigh(du,u,p,t)
    T, M, m, L, num_bath_particles = p
    Utag(x) = -T * (L - 2* x) / (x*(L-x)) 
    Gamma(x) = num_bath_particles * (1/x + 1/(L-x)) * sqrt(8*π * m * T)
    du[1] = u[2]
    du[2] = (-Utag(u[1]) - Gamma(u[1]) * u[2]) / M
end

function σ_rayleigh(du,u,p,t)
    T, M, m, L, num_bath_particles = p
    Gamma(x) = num_bath_particles * (1/x + 1/(L-x)) * sqrt(8*π * m * T)
    du[1] = 0
    du[2] = sqrt(2* abs(Gamma(u[1])) * T) / M
end

T = 1.
M = 10^6
m = 1.
L = 10.
num_bath_particles = 1
p = [T, M, m, L, num_bath_particles]
x0 = 2.5
v0 = 0.0
prob_sde_rayleigh = SDEProblem(rayleigh, σ_rayleigh, [x0, v0], (0.0,10^6), p)
ensembleprob = EnsembleProblem(prob_sde_rayleigh)

dt = 1000.
sim = solve(ensembleprob, EnsembleThreads(), trajectories=1000, saveat=dt) 

using DifferentialEquations.EnsembleAnalysis
summ = EnsembleSummary(sim, 0:dt:10^6)
plot(summ,labels="Middle 95%", idxs=1)

and for the 2nd plot is

dt = 1000
t_to_run = 10^6
t_arr = 1:dt:t_to_run
m_series, v_series = timeseries_point_meanvar(sim, 1:dt:t_to_run)
x_mean= zeros(length(dt:dt:t_to_run))
x_var= zeros(length(dt:dt:t_to_run))
v_mean= zeros(length(dt:dt:t_to_run))
v_var = zeros(length(dt:dt:t_to_run))
for i=1:length(dt:dt:t_to_run)
    x_var[i] = v_series.u[i][1]
    x_mean[i] = m_series.u[i][1]
    v_var[i] = v_series.u[i][2]
    v_mean[i] = m_series.u[i][2]
end
plot(x_mean)

Why is there such inconsistency?

2

I think the top plot is displaying the median and the bottom plot is displaying the arithmetic mean. For multiplicative noise processes these tend to diverge.

cmp

4 Likes
3

Yes, IIRC the default for the summary plot is to use quantiles, so it’ll be showing the median and the 95% quantiles, which in some cases of course can be very different from the mean.

4

@josuagrw thanks very much. I did not realize the summary gives the median.
@ChrisRackauckas what is IIRC?

1 Like
5

IIRC — if I recall correctly

1 Like
6

Just as a sidenote: If your model is correct, this could be a quite interesting scientific result. It would mean that there are very rare cases which have such intensity that they overshadow the more frequent “normal” cases. Nassim Taleb writes extensively on the risk implications for such systems.

1 Like
7

Yes, that’s a good point. Since you told me that I was actually looking at the median, I am trying to figure out if I have a problem with my model.

The coordinate I am looking at should be bounded in some range (in this case between 0 and 10), and this is modeled by an infinite force at x=0, 10 which rises in a continuous though fast manner. I suspect that the solver manages to jump over this barrier and then go wild.
I am looking at very large times, so the solver might take too coarse timesteps, even though I was thinking that it should find appropriate timesteps.

8

You might need Stratonovich integration to describe this process correctly…

1 Like
9

thanks for the tip! I will try to figure it out.

10

I don’t think this is the case. For any Itô SDE, one can write down an equivalent Stratonovich SDE, and similarly for the reverse (perhaps assuming mild technical hypotheses). For a brief reference, see here.

Maybe the modeling would be more intuitive/practical in the Stratonovich interpretation, however.

3 Likes
11

What I mean is, there might be a mismatch between the form of the model (might be stratonovich, common in physics) and the numerical solver (which is ito by default).

12

My apologies, I misunderstood your comment!

1 Like
13

I tried some of the Stratonovich solvers, like SRA and SOSRA etc., (without putting too much thought into it yet) and some did yield a better result. This means that for some ensembles the mean explodes and some do not. I will further investigate this.

14

That’s interesting. I want to stress @patrick s point: Switching between Ito and Stratonovich integration changes the meaning of the parameters in your model. So you need to make sure beforehand in which form your SDE is formulated and then choose the correct type of solver for that formulation.

The article @patrick linked to is indeed extremely enlightening on that topic.

1 Like
15

SRA and SOSRA aren’t not for Ito or Stratonovich necessarily, just additive noise, and additive noise is the same in both interpretations.

2 Likes