Hi, I’ve been working with this system of differential equations
function JJ!(du,u,p,t)
I,g=p
du[1] = - gu[1] - 2u[3] + I
du[2] = -g*u[2] + u[1]u[3]
du[3] = -gu[3] +(1 -u[2])*u[1]
end
I solve it using the DifferentialEquations library specifically I use:
function get_solutions(omega,g,I,x0,y0,T,method=AutoTsit5(Rosenbrock23()))
tspan = (0.0,T)
prob = ODEProblem(JJ!,[omega,x0,y0],tspan,[I,g])
x=1e-5
sol=solve(prob, method, abstol = x, reltol = x,)
return sol
end
As I change x which is my error the solutions change dramatically in the long time behaviour, If I just use the default I get what I expect but if I increase it, the behaviour I expect (oscillations) are no longer there, also If I keep default and use some other method such as Vern7 then it also changes, what’s a good guide to tell if the long time behaviour is real or a numerical artifact?
Convergence. What happens as the tolerance goes to zero?
I get the same behavior just highly attenuated. If I set up atol=rtol=1e-17 then I get oscillations of around 1e-14 for AutoTsit5(Rosenbrock23()) and 1e-15 for Vern7(), while for atol=rtol=1e-10 I get oscillations of about 1e-8 for Autosit and 1-e11 for Vern7(). The shape of the oscillations is also different. I think the oscillations might be there because they seem to be some orders of magnitude below the tolerance when it goes to zero but at the same time, it’s like their amplitude is going to zero so I’m kind of confused. Should I increase maxiters, to see what happens at atol=rtol=1e-30 for example?
What about with Rodas4, CVODE_BDF, or radau? Higher order methods for stiff ODEs is usually a good test for correctness.
Sorry for the delay in my response, Using Rodas4 and RaudauIIA5 I get different things in long times I’m setting atol=rtol=1e-15
This is with interpolation, using
t0 = 4000 # Start time
fs = 400
tmax = 5000 # End time
t = t0:1/fs:tmax;
A = sol(t)
I say that because I feel interpolation might be part of the problem and I don’t know how to turn it off. for CVODE_BDF I get
I couldn’t get lower atol=retol because I get [CVODES ERROR] CVode
At t = 0, too much accuracy requested. my coal is to verify that oscillations are there in the long time limit
If atol=rtol is 1e-5 for example, then I get no oscillations for Rodas or radau but get them for CVODE_BDF, should I trust the CVODE_BDF and think the oscillations are really there? or is there any further test I should try to do
