Warning dt<dtmin in ODE Solver not solved with FAQ

General Usage
1

Dear all, I have run into the following issue, which is described in details also here: Frequently Asked Questions · DifferentialEquations.jl (sciml.ai).

I have tried all suggestions on the link and made a great amount of small modifications to test alternatives. Also, the same exact implementation presents no instability in MATHEMATICA, so I strongly believe the model is correct.

Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase C:\Users\18143\.julia\packages\SciMLBase\L7Nun\src\integrator_interface.jl:345

Here is the main part of my code relevant to this:

## 
##Define the parameters of the dynamics
dt_step=0.10
time=collect(dt_step:dt_step:4.0)
initialcond=collect(0.1:0.2:0.9)
N=4; M=4; T=4; L=1; Nsamples=(N+1)*(M+1)

##Problem-related parameters 
β=0.5; ξ=0.5; α=-2; C=2;

### --- Functions --- ###
function uₜₑᵣₘₛ(i,j,t,x₀)
    uₜₑᵣₘₛ = aₘₙ[i,j]*cos((i-1)*π*x₀/L)*cos((j-1)*π*t/T)
end

function uₐₚₚ(t,x₀)
    uₐₚₚ = sum(uₜₑᵣₘₛ(i,j,t,x₀) for i=1:(M+1), j=1:(N+1))
end

### --- More Functions --- ###
using DifferentialEquations

function ϕᵣₑₐₗ(x₀)
    f(ϕ,p,t) = β*ϕ*(1-ϕ)*(uₐₚₚ(t,x₀)-ξ)
    ϕ₀=x₀
    tspan=(0.0,4.0)
    prob=ODEProblem(f,ϕ₀,tspan)
    sol=solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) #Tested many of these already
    return sol
end

#Numerical Integration of Functional - Trapezoidal Rule
function quad_trap(f,a,b,NN,x₀) 
    h = (b-a)/NN
    int = h * ( f(a,x₀) + f(b,x₀) ) / 2
    for k=1:NN-1
        xk = (b-a) * k/NN + a
        int = int + h*f(xk,x₀)
    end
    return int
end

function RealFunctional(t,x₀)
    α*ϕᵣₑₐₗ(x₀)(t)*uₐₚₚ(t,x₀) + C*uₐₚₚ(t,x₀)^2
end

function RealFunctionalCost(x₀)
    Int = quad_trap(RealFunctional,0,T,size(time)[1],x₀)
    return Int
end

function RealJₐ()
       return sum(RealFunctionalCost(x₀) for x₀=initialcond)
end

#Gradient Descent 
aₘₙ=rand(Uniform(-1,1),(M+1),(N+1)); ϵ=10e-2;
ApproxRealGrad=Matrix{Float64}(undef,(M+1),(N+1))

#Gradient Computation -- I know this is lousy but it needs to be like this
    FunctionalI=RealJₐ()
    println(FunctionalI)
    a_δ=aₘₙ
        for i=1:(M+1), j=1:(N+1)
            aₘₙ=a_δ
            aₘₙ[i,j] = a_δ[i,j] + ϵ
            FunctionalF=RealJₐ()
            ApproxRealGrad[i,j] = (FunctionalF-FunctionalI)/ϵ
        end    
    aₘₙ=a_δ

τ=0.1
    #Descent 
    aₘₙ -= τ*ApproxRealGrad

There is a loop in the above omitted here, because I am running iteration by iteration “by hand” in Jupyter.

Also, I know there would be many alternatives to the computation of this derivative, but for our application this is a necessary step. The same step in the descent step worked fine in MATHEMATICA as well.

Please help me, mainly with the ODE ϕᵣₑₐₗ(x₀) , which is the source of this.

I spent a lot of time already trying to get out of this warning. Thank you a lot for your time and help.

Gabriel

2

Following:

Was the instability with all solvers? Including CVODE and the Fortran methods? Those are the same solvers you can use from Mathematica, so you should get the same results if you choose them.

1 Like
3

Thank you again for your support. I made it work with the Sundials.jl. Now I have exactly the same output as in Mathematica.

Gabriel

4

Could you post an example that can be explored? I’m curious why for example FBDF() wouldn’t work if Sundials does. That’s worth fixing if that’s the case.