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

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
``````

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?

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)
``````

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
``````

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
``````

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

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.

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

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

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);
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))
``````

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.

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