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!