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

**URL:** <https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, sde\
**Created:** [October 31, 2020, 10:52pm UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392 "2020-10-31T22:52:41Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [October 31, 2020, 10:52pm UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/1 "2020-10-31T22:52:41Z")

</div>

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

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

```

for certain inputs. See this example:

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

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

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

```

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

 ![sde_sim_problem](https://global.discourse-cdn.com/julialang/original/3X/a/5/a552725be503dd6471f7d1b47ebd590a9fefd8ea.png)

This is the output of `Pkg.status()`

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 1, 2020, 12:19am UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/2 "2020-11-01T00:19:27Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [November 1, 2020, 12:26am UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/3 "2020-11-01T00:26:55Z")

</div>

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

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

```

![dt_over_t](https://global.discourse-cdn.com/julialang/original/3X/8/2/826e59b0fc73816c284ed1e4dc00a79e93728035.png)

The last 20 values of `sol.t`

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

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

```

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [November 1, 2020, 12:31am UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/4 "2020-11-01T00:31:19Z")

</div>

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:

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

```

gives

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

```

and

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

```

gives

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 1, 2020, 1:36am UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/5 "2020-11-01T01:36:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [November 1, 2020, 4:03pm UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/6 "2020-11-01T16:03:25Z")

</div>

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

I tried printing

```julia
sol.destats

```

which gives the output

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

 ![maxiterserror](https://global.discourse-cdn.com/julialang/original/3X/7/d/7d5fb80f06ceb97db902a76e78c8c4672f91be48.png)  
I also did the dt length comparison of the two solutions.  
 ![dt_over_t](https://global.discourse-cdn.com/julialang/original/3X/c/d/cd5c5342eb29c5dc3c688c8e086644bf08d4a920.png)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 1, 2020, 5:48pm UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/7 "2020-11-01T17:48:42Z")

</div>

> [@Torkel](#):
>
> In what field of the integrator can I find that information?

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

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [November 1, 2020, 6:06pm UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/8 "2020-11-01T18:06:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [November 2, 2020, 5:57pm UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/9 "2020-11-02T17:57:48Z")

</div>

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.

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

```

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

```

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

```

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

```

 ![trajs](https://global.discourse-cdn.com/julialang/original/3X/0/f/0ff963041b266d21121111d075ba41372712f14a.png)

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

```

 ![dts](https://global.discourse-cdn.com/julialang/original/3X/4/f/4f762034d87962db467b969a0039d18f5368f409.png)

Happy to provide any more info!

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [November 14, 2020, 3:44pm UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/10 "2020-11-14T15:44:17Z")

</div>

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

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

 ![sol_fails](https://global.discourse-cdn.com/julialang/original/3X/0/6/06e15e237256eee3eea7c35e7f4bbf76d15422fa.png)

dts:

 ![dt_fails](https://global.discourse-cdn.com/julialang/original/3X/e/8/e8b77cdbe003cbe2366c5c858415bf57332f3781.png)

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`

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

```

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

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [November 14, 2020, 3:50pm UTC](https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392/11 "2020-11-14T15:50:06Z")

</div>

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