Problem for a system of stochastic differential equations (SDE) with multiple Wiener processes

[canturk]

I attempted to solve a system of SDE as in the form below. But I am struggeling with a persisting error message as shown below.

A replicated SDE is as follows:

Here A, B, C, and D are sparse matrices with size NxN; u is a vector of length N; and dw[i] is ith Wiener process forall i in {B,C,D}.

For implementation purpose, SDE will be restructured as in the form

The following is an example code (based on a post that generates a WARNING "There is either an error in your model specification or the true solution is unstable".

using DifferentialEquations, SparseArrays
# Parameters
N = 10; 
tspan = (0.0, 10.0);
ntraj = 2;
# Matrices and initial value
A = sprand(N, N, 0.5); B = sprand(N, N, 0.5); C = sprand(N, N, 0.5); D = sprand(N, N, 0.5)
u0= rand(N);
# Deterministic and stochastic functions
function f!(Au, u, p, t) 
    Au .= A*u
    nothing
end
function g!(Gu, u, p, t)
    Gu[:,1] = B * u
    Gu[:,2] = C * u
    Gu[:,3] = D * u
    nothing
end;
# Construct SDE problem with three Wiener processes
prob = SDEProblem(f!, g!, u0, tspan, noise_rate_prototype=sparse(zeros(N,3)));
# Solve the equation of motion for an arbitrary number of trajectories
# sol = solve(EnsembleProblem(prob), LambaEM(), save_everystep=false, EnsembleThreads(),trajectories=ntraj)
sol = solve(EnsembleProblem(prob), LambaEulerHeun(), save_everystep=false, EnsembleThreads(),trajectories=ntraj)

Here is the warning message:

I was wondering why this error message pops up: Is it because of an error in the model or anyting else?

Any help or suggestion is appreciated.

[contradict]

I think it may actually be unstable. I tried generating several matrices with the method you described and always get eigenvalues > 1.0, so I expect the solution to grow. Try multiplying all your matrices by 1.0/N.

julia> using LinearAlgebra, SparseArrays, Arpack

julia> N = 10
10

julia> A = sprand(N,N, 0.5);

julia> λ, ϕ = eigs(A);

julia> abs.(λ)
6-element Vector{Float64}:
 2.2336327664915077
 1.1642719859775836
 1.1642719859775836
 0.5622555935660684
 0.5622555935660684
 0.45944566554197847

[BLI]

The eigenvalues should be negative in continuous time.

[contradict]

Oops, you are correct! I should have been looking at

ulia> real.(λ)
6-element Vector{Float64}:
  2.3685798031036582
  1.2284625269151874
 -0.3457237484386846
 -0.3457237484386846
 -0.7113750263386742
 -0.7113750263386742

Which gives the same answer.

[BLI]

…ah, I should have been more precise: negative real part, assuming linear deterministisk part and no state dependence on stochastic part…

[canturk]

Thanks for the suggestion but multiplying every matrix with 1/N didn’t help for a convergent solution. However I think if we would be able to normalize the solution vector u on very step during the evolution might fix the issues. So I was wondering whether it is possible to force the SDE solver to normalize in every step.