Strange MaxIters (numeric instability?) occured when solving certain SDE

Gaussia · October 31, 2020, 10:52pm

I’ve been solving this equation for a lot of different parameters, but suddenly I am getting occasional

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase /home/torkelloman/.julia/packages/DiffEqBase/V7P18/src/integrator_interface.jl:329

for certain inputs. See this example:

using StochasticDiffEq

function sys_f(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1] = v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1) - σ
    du[2] = τ*(σ-A)
end
function sys_g(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1,1] = η*sqrt(max(v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1),0.))
    du[1,2] = -η*sqrt(max(σ,0.))
    du[1,3] = 0
    du[1,4] = 0
    du[2,1] = 0
    du[2,2] = 0
    du[2,3] = η*sqrt(max(τ*σ,0.))
    du[2,4] = -η*sqrt(max(τ*A,0.))
end


u0 = [0.025, 0.025]
tspan = (0.,200.)
p = [28.48035868435802, 0.2527461586781728, 0.05, 0.025, 3.0, 0.1]

sde_prob = SDEProblem(sys_f,sys_g,u0,tspan,p,noise_rate_prototype=zeros(2,4));

sol = solve(sde_prob,ImplicitEM(),seed=4963697075065472807);

Solving without this seed, things are mostly fine:

for repeat = 1:1000
    sol = solve(sde_prob,ImplicitEM())
end

yields maybe only 1 such error per run.

Even increasing maxiters a bit, I receive no benefit:

sol = solve(sde_prob,ImplicitEM(),seed=4963697075065472807,maxiters=10000000);

still errors. Finally, I am attaching an image of bad/good simulation.

This is the output of Pkg.status()

  [717857b8] DSP v0.6.8
  [0c46a032] DifferentialEquations v6.15.0
  [5ab0869b] KernelDensity v0.6.1
  [91a5bcdd] Plots v1.6.8
  [f27b6e38] Polynomials v1.1.10
  [10745b16] Statistics

ChrisRackauckas · November 1, 2020, 12:19am

That plot of the “bad” one doesn’t show much. Why did it hit max iters? Smaller dt? Why is dt smaller? Looks like the relative noise would increase?

Gaussia · November 1, 2020, 12:26am

Seems dt is getting small, this is a plot of:

sol = solve(sde_prob,ImplicitEM(),seed=4963697075065472807);
dts = diff(sol.t)
plot(sol.t[2:end],dts,framestyle=:box)

dt_over_t

The last 20 values of sol.t

3.0696546597633163
 3.0696553650908487
 3.0696561585843227
 3.0696568640351365
 3.069657657667302
 3.069658363134364
 3.069659156784809
 3.0696598622394653
 3.069660655875954
 3.0696613614302923
 3.069662155178923
 3.0696628609243257
 3.0696636548879033
 3.0696643602892197
 3.0696651538657007
 3.0696658593007053
 3.0696666529150853
 3.0696673582528358
 3.069668151757805
 3.069668856962127
 3.0696696503169894

the last 20 values of dts

7.934824788335959e-7
 7.053275323798402e-7
 7.934934740383426e-7
 7.054508137649407e-7
 7.936321653190248e-7
 7.054670621009507e-7
 7.936504449190807e-7
 7.054546564688735e-7
 7.936364885274827e-7
 7.055543385092733e-7
 7.937486308229325e-7
 7.057454025627408e-7
 7.939635775500165e-7
 7.054013164697892e-7
 7.935764809730017e-7
 7.054350046331592e-7
 7.936143799902595e-7
 7.053377504284697e-7
 7.935049692875396e-7
 7.052043220490134e-7
 7.933548622496289e-7

Gaussia · November 1, 2020, 12:31am

Noise does not seem to dominate either. I am a bit unsure about this, but basically checking the values of f and g at the final point of the solution:

u_in = zeros(2)
sys_f(u_in,sol.u[end],p,0.)
u_in

gives

2-element Array{Float64,1}:
 0.5475974457463257
 0.023042909797453985

and

u_in = zeros(2,4)
sys_g(u_in,sol.u[end],p,0.)
u_in

gives

2×4 Array{Float64,2}:
 0.101223  -0.0690655  0.0         0.0
 0.0        0.0        0.0154435  -0.00284127

ChrisRackauckas · November 1, 2020, 1:36am

Can you poke into the integrator and check which thing is causing rejection? Is EEst > 1 or is it the Newton iterations diverging?

Gaussia · November 1, 2020, 4:03pm

In what field of the integrator can I find that information?

I tried printing

sol.destats

which gives the output

DiffEqBase.DEStats
Number of function 1 evaluations:                  2450
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    20862
Number of linear solves:                           0
Number of Jacobians created:                       2450
Number of nonlinear solver iterations:             0
Number of nonlinear solver convergence failures:   0
Number of rootfind condition calls:                0
Number of accepted steps:                          0
Number of rejected steps:                          0

I also tried printing both solutions, but over the same interval t=(0.,3.). The faulty one looks weird throughout.

I also did the dt length comparison of the two solutions.

ChrisRackauckas · November 1, 2020, 5:48pm

It’s not in the integrator. You need to @show in the source for that.

Gaussia · November 1, 2020, 6:06pm

Ok, still not really sure how to to “@show in the source”. What is the source? The actual sys_f function?

Gaussia · November 2, 2020, 5:57pm

I got access to a rather large number of parameters for which this happens. Most notably, the last parameter scales the noise. This error happens for this parameter is 0.001 (very small). Here the noise is essentially negligible.

using StochasticDiffEq

function sys_f(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1] = v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1) - σ
    du[2] = τ*(σ-A)
end
function sys_g(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1,1] = η*sqrt(max(v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1),0.))
    du[1,2] = -η*sqrt(max(σ,0.))
    du[1,3] = 0
    du[1,4] = 0
    du[2,1] = 0
    du[2,2] = 0
    du[2,3] = η*sqrt(max(τ*σ,0.))
    du[2,4] = -η*sqrt(max(τ*A,0.))
end

u0 = [0.025, 0.025]
tspan = (0.,200.)
p = [6.7341506577508214, 1.2224984640507832, 0.05, 0.025, 3.0, 0.001]

sde_prob = SDEProblem(sys_f,sys_g,u0,tspan,p,noise_rate_prototype=zeros(2,4));

sol1 = solve(sde_prob,ImplicitEM(),seed=5105556673070507935); sol = sol1;
sol2 = solve(sde_prob,ImplicitEM());

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase /home/torkelloman/.julia/packages/DiffEqBase/V7P18/src/integrator_interface.jl:329

using Plots
gr(); default(fmt = :png);
p1 = plot(sol1,framestyle=:box,title="Bad stuff",xlimit=(0.,sol1.t[end]))
p2 = plot(sol2,framestyle=:box,title="Everything is fine",xlimit=(0.,sol1.t[end]))
p3 = plot(sol2,framestyle=:box,title="Everything is fine",xlimit=(0.,sol2.t[end]))
plot_trajectories = plot(p1,p2,p3,size=(1200,400),layout=(1,3))

p1 = plot(sol1.t[2:end],diff(sol1.t),framestyle=:box,xlimit=(0.,sol1.t[end]),title="Bad stuff")
p2 = plot(sol2.t[2:end],diff(sol2.t),framestyle=:box,xlimit=(0.,sol1.t[end]),title="Everything is fine")
p3 = plot(sol2.t[2:end],diff(sol2.t),framestyle=:box,xlimit=(0.,sol2.t[end]),title="Everything is fine")
plot_dts = plot(p1,p2,p3,size=(1200,400),layout=(1,3))

Happy to provide any more info!

Gaussia · November 14, 2020, 3:44pm

I have tried investigating this over several solvers. Looking at the initial problem I try:

sol_ImplicitEM = solve(sde_prob,ImplicitEM(),seed=4963697075065472807);
sol_STrapezoid = solve(sde_prob,STrapezoid(),seed=4963697075065472807);
sol_SImplicitMidpoint = solve(sde_prob,SImplicitMidpoint(),seed=4963697075065472807);
sol_ImplicitEulerHeun = solve(sde_prob,ImplicitEulerHeun(),seed=4963697075065472807);
sol_ISSEM = solve(sde_prob,ISSEM(),seed=4963697075065472807);
sol_ISSEulerHeun = solve(sde_prob,ISSEulerHeun(),seed=4963697075065472807);

Of these, sol_ImplicitEM, sol_STrapezoid , sol_SImplicitMidpoint all fails, while
ImplicitEulerHeun, ISSEM, ISSEulerHeun succeeds. See plots (know it is not super helpful, but best I got to show how they fail.
Solutions:

dts:

Finally, I try scanning for other seeds which causes failure in the remaining three. Did not find any for ImplicitEulerHeun or ISSEulerHeun, however, ISSEM fails for seed=15088718199444158177

sol_ISSEM = solve(sde_prob,ISSEM(),seed=15088718199444158177);

Looking at this seed specifically, ImplicitEM and SImplicitMidpoint also fails, while the remaining three succeeds.

Gaussia · November 14, 2020, 3:50pm

Trying to find a pattern, it seems like I can find errors for the Ito methods, but not for the Stratonovich methods.