Skip to content

New release breaks initialization of simple ODE #920

@hersle

Description

@hersle

With the upgrade from 2.71.4 to 2.72.0 (in particular commit 204987b), this example stops working:

using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables a(t)
@named sys = ODESystem([
    D(a) ~ (1e-5/a^2) * a^1
], t; initialization_eqs = [
    a ~ D(a) * t
], guesses = [
    a => 1
])
sys = structural_simplify(sys)
prob = ODEProblem(sys, [], (1e-5, 1.0), [])
sol = solve(prob, Tsit5())
@test sol[a][1]  1e-5^(3/2) # anal solution of D(a) = √(1e-5) --> a = √(1e-5)*t + C

The solve(prob, Tsit5()) fails with InitialFailure:

┌ Warning: Automatic dt set the starting dt as NaN, causing instability. Exiting.
└ @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/Ifi3S/src/solve.jl:649
retcode: InitialFailure
Interpolation: specialized 4th order "free" interpolation
t: 1-element Vector{Float64}:
 1.0e-5
u: 1-element Vector{Vector{Float64}}:
 [0.0]

However, solve(prob.f.initialization_data.initializeprob) succeeds with

retcode: Success
u: 1-element Vector{Float64}:
 3.162277661949986e-8

Metadata

Metadata

Labels

bugSomething isn't working

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions