Inconsistency when solving system of SDE with diagonal noise

Hello,
I noticed a strange behaviour of Julia when trying to solve the following system of SDE:

\begin{align} du_1 &= \left(\frac{1}{1+c^2}\right) \left( -\frac{u_4 - 2u_3+ u_4}{2d^2} + (u_1^2 + u_3^2)u_3 - c \left( -\frac{u_2 - 2u_1 + u_2}{2d^2} + (u_1^2 + u_3^2)u_1 \right) \right)dt +dW_1\\ du_2 &= \left(\frac{1}{1+c^2}\right) \left( -\frac{u_3 - 2u_4+ u_3}{2d^2} + (u_2^2 + u_4^2)u_4 - c \left( -\frac{u_1 - 2u_2 + u_1}{2d^2} + (u_2^2 + u_4^2)u_2 \right)\right)dt+dW_2 \\ du_3 &= \left(-\frac{1}{1+c^2}\right) \left( -\frac{u_2 - 2u_1+ u_2}{2d^2} + (u_1^2 + u_3^2)u_1 + c \left( -\frac{u_4 - 2u_3 + u_4}{2d^2} + (u_1^2 + u_3^2)u_3 \right) \right)dt+dW_3 \\ du_4 &= \left(-\frac{1}{1+c^2}\right) \left( -\frac{u_1 - 2u_2+ u_1}{2d^2} + (u_2^2 + u_4^2)u_2 + c \left( -\frac{u_3 - 2u_4 + u_3}{2d^2} + (u_2^2 + u_4^2)u_4 \right) \right)dt+dW_4 \\ \end{align}

where (dW_1,dW_2,dW_3,dW_4) are independent Winer increments. So this falls under the category of system of SDE with diagonal noise.
Now, the way I coded it is the following:

using DifferentialEquations

#parameters 
c = 0.03          
d= 0.3         
p=(c=c,d=d)

#deterministic function F
function F!(du,u,p,t)
        du[1]=(1/(1+p.c^2)) * ( -(u[4]-2*u[3]+u[4])/(2*p.d^2)
                             +(u[1]^2+u[3]^2)*u[3]
                             - p.c*(-(u[2]-2*u[1]+u[2])/(2*p.d^2) +(u[1]^2+u[3]^2)*u[1]))
        du[2]=(1/(1+p.c^2)) * ( -(u[3]-2*u[4]+u[3])/(2*p.d^2)
                                                  +(u[2]^2+u[4]^2)*u[4]
                                                  -p.c*(-(u[1]-2*u[2]+u[1])/(2*p.d^2) +(u[2]^2+u[4]^2)*u[2]))
        du[3]=(-1/(1+p.c^2)) * ( -(u[2]-2*u[1]+u[2])/(2*p.d^2)
                                                  +(u[1]^2+u[3]^2)*u[1] 
                                                  +p.c*(-(u[4]-2*u[3]+u[4])/(2*p.d^2) +(u[1]^2+u[3]^2)*u[3]))                                          
			
        du[4]=(-1/(1+p.c^2)) * ( -(u[1]-2*u[2]+u[1])/(2*p.d^2)
			                              +(u[2]^2+u[4]^2)*u[2] 
			                              +p.c*(-(u[3]-2*u[4]+u[3])/(2*p.d^2) +(u[2]^2+u[4]^2)*u[4]))                                              
		
end	

#noise function 
function G!(du, u, p, t)        	
	for i in 1:4
		du[i]=1.0
	end		
end

#solution
u_0=zeros(4)
tspan=(0.0,5.0)
problem=SDEProblem(F!,G!,u_0,tspan,p)
sol=solve(problem,save_everystep=false)

This seems to work, as I get the following message:

retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 5.0
u: 2-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 0.0]
 [-0.2810927281000788, -0.3205507707167091, -0.8684132010707524, -1.131861129900748]

However now comes the strange thing. Let’s rewrite the noise term in the following way (which I believe should be equivalent to the previous!)

function G!(du, u, p, t)
 	for i in 1:4
        for j in 1:4
            if (i==j)
                du[i,j]=1.0
            else
                du[i,j]=0.0
            end
        end
    end              	
			
end

#solution
u_0=zeros(4)
tspan=(0.0,5.0)
problem=SDEProblem(F!,G!,u_0,tspan,p,noise_rate_prototype=zeros(4,4))
sol_r=solve(problem,save_everystep=false)

So basically what I changed here is that now G is the 4\times4 identity matrix, which I specified to the solver with:

noise_rate_prototype=zeros(4,4)

and dW is a vector containing 4 independent noises (dW_1,dW_2,dW_3,dW_4), such that the matrix-vector product G\cdot dW correctly gives the change of the solution due to the stocastic part of the equation. (this is also explained in the documentation https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/ , in particular example 4).
The problem with this is that the code with this modification gives the following instability warning:

┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=2.3460875073396528, and step error estimate = 0.754176907452515. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/Dwomw/src/integrator_interface.jl:619
retcode: DtLessThanMin
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 2.3460875073396528
u: 2-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 0.0]
 [-3.6756952266070644, -1.7968454884186287e6, -0.5435519681145712, 1.19293585662674e7]

However, as I said, I thought that this second way was equivalent to the first, which didn’t give any warning. So what is going on? I ask this because in the end I’m interested in implementing non-diagonal noise using this second approach (with a G matrix non diagonal), but already for G diagonal I run into this instability problem and this inconsistency with the previous approach.
Thanks!

The higher order methods for diagonal noise are much better than those for non-diagonal noise. If you use non-diagonal noise, it needs to switch to an algorithm that has lower stability, worse error estimators, and worse performance. If you have a diagonal noise SDE you really want to make sure you define it as diagonal to get all of the benefits of the much better solvers: non-diagonal SDEs are mathematically a hundred times harder than diagonal SDEs!

1 Like