I am using SciML to try to integrate a stiff model with small noise. It is giving me tons of grief though, and I am not sure how to approach it. I even tried setting tolerances to very small values & increasing maxiter to several orders of magnitude above the default.
As a test, I attempted to integrate it using zero noise, since I know the deterministic system can be integrated stably. However, as soon as I switch from the deterministic system with AutoVern7(Rodas5()) to the system with zero noise amplitude using ImplicitRKMil() or SKenCarp() it interrupts the integration.
With ImplicitRKMil() I get an error like: Interrupted. Larger maxiters is needed.
With SKenCarp(): At t=1843.1833699767562, dt was forced below floating point epsilon 1.1758922554404307e-13, and step error estimate = 0.20233022763024602. Aborting. There is either an error in your model specification or the true solution is unstable.
Here is an MWE:
using DifferentialEquations
using StochasticDelayDiffEq
using Plots
const w=[-0.0974439293080174, -0.0972907033232282, 0.05008580780639362, 31.382671443379348, 1.1925811288818222, 45.142984062837634]
function eqns!(du,u,h,p,t)
global w
eps1=0.09
eps2=0.003
FREQ=0.9
delta=p[1]
CNORM=2.0*pi/FREQ
u1 =u[1]
u2 = u[2]
du[1]=CNORM*(u2 - delta*cos(2*pi*t + pi/2))
du[2]=CNORM*(1.0 - eps1*u2 - exp(u1) - eps2*exp(u1)*u2)
du[1]+=w[1]*h(p,t-w[4])[1]+w[2]*h(p,t-w[5])[1]+w[3]*h(p,t-w[6])[1]
du[2]+=w[1]*h(p,t-w[4])[2]+w[2]*h(p,t-w[5])[2]+w[3]*h(p,t-w[6])[2]
end
function g!(du,u,h,p,t)
du[1]=p[2]*cos.(2.0*pi*t)
du[2]=0
end
#this works fine.
h_=function(p,t) [-5.7734959978E+01,4.6557062215E-01] end
stable_problem = DDEProblem(eqns!, h_(0,0),h_, (0.0,5000), (6.08,0);
constant_lags = (w[4],w[5],w[6]),dependent_lags=[]
)
alg=MethodOfSteps(AutoVern7(Rodas5()))
sol=solve(stable_problem,alg,dense=true)
png(plot(sol),"stable.png") #these results are sensible
#this fails at time t=1843.1833... with a dt forced below floating point epsilon warning
#note that the noise amplitude is zero
#if I increase the noise slightly, it is the same thing.
same_problem = SDDEProblem(eqns!, g!, h_(0,0), h_, (0.0,5000), (6.08,0);
constant_lags = (w[4],w[5],w[6]),dependent_lags=[])
#this tells me I need larger maxiters and then fails also
new_alg=SKenCarp()
sol2=solve(same_problem,new_alg,dense=false,saveat=0.01)
png(plot(sol2), "unstable.png") #if I decrease the integration time to 1000 this will produce some sort of nonempty figure, but it's way off
another_alg=ImplicitRKMil()
sol3=solve(same_problem,another_alg,dense=false,saveat=0.01)
png(plot(sol3), "unstable2.png")
Any ideas?