DifferentialEquations.jl
DifferentialEquations.jl copied to clipboard
Published
20 hours ago •
SciML
SciML
"ODEProblem: Warming:Instability detected. Aborting" when using Mutable Arrays with ManifoldProjection
This
using DifferentialEquations, StaticArrays
function harmonic!(du,u,p,t)
du[1] = u[2]
du[2] = -u[1]
end
u0 = @MArray [1.0, 0.0]
function constantenergymanifold(resid,u,p,t)
resid[1] = u[2]^2 + u[1]^2 - 1.0
resid[2] = 0.0
end
constantenergy_cb = ManifoldProjection(constantenergymanifold)
prob = ODEProblem(harmonic!, u0, (0., 100.), nothing, callback=constantenergy_cb)
sol = solve(prob)
gives warning.
┌ Warning: Instability detected. Aborting
└ @ DiffEqBase C:\Users\DeadScholar\.julia\packages\DiffEqBase\8fdjX\src\integrator_interface.jl:339
retcode: Unstable
Interpolation: Automatic order switching interpolation
t: 27-element Array{Float64,1}:
0.0
0.0009990009990009992
0.0009990009990009992
0.010989010989010992
0.010989010989010992
0.07985922249873038
0.07985922249873038
0.2403882280626971
0.2403882280626971
0.48125583199780264
0.48125583199780264
0.7872724750204353
0.7872724750204353
⋮
1.6302664787742993
1.6302664787742993
2.1040017944316363
2.1040017944316363
2.6513320154468447
2.6513320154468447
3.2715753683812103
3.2715753683812103
3.8874031049079534
3.8874031049079534
4.57110644570036
4.57110644570036
u: 27-element Array{MArray{Tuple{2},Float64,1,2},1}:
[1.0, 0.0]
[0.9999995009985435, -0.000999000832833342]
[0.9999995009985435, -0.000999000832833342]
[0.9999396214263468, -0.010988789821184267]
[0.9999396214263468, -0.010988789821184267]
[0.9968129466114029, -0.07977436592568171]
[0.9968129466114029, -0.07977436592568171]
[0.9712456180254766, -0.23807970964822986]
[0.9712456180254766, -0.23807970964822986]
[0.8864142948207253, -0.4628927349710369]
[0.8864142948207253, -0.4628927349710369]
[0.7057801319307301, -0.7084309070362929]
[0.7057801319307301, -0.7084309070362929]
⋮
[-0.05943636038391344, -0.9982319821376294]
[-0.05943636038391344, -0.9982319821376294]
[-0.508298495902857, -0.8611808133891722]
[-0.508298495902857, -0.8611808133891722]
[-0.8822130263524072, -0.4708509609744615]
[-0.8822130263524072, -0.4708509609744615]
[-0.9915642536268217, 0.12963020212303705]
[-0.9915642536268217, 0.12963020212303705]
[-0.7345269620571533, 0.6785843850919284]
[-0.7345269620571533, 0.6785843850919284]
[-0.14078005159723983, 0.9900486388892646]
[NaN, NaN]
Changing u0 = @MArray [1.0, 0.0] to u0 = [1.0, 0.0] makes it work
retcode: Success
Interpolation: Automatic order switching interpolation
t: 253-element Array{Float64,1}:
0.0
0.0009990009990009992
0.0009990009990009992
0.010989010989010992
0.010989010989010992
0.07985922249873038
0.07985922249873038
0.2403882280626971
0.2403882280626971
0.48125583199780264
0.48125583199780264
0.7872724750204353
0.7872724750204353
⋮
95.74825196316992
95.74825196316992
96.63761832689438
96.63761832689438
97.46808053105032
97.46808053105032
98.34371732963653
98.34371732963653
99.19194600937247
99.19194600937247
100.0
100.0
u: 253-element Array{Array{Float64,1},1}:
[1.0, 0.0]
[0.9999995009985435, -0.000999000832833342]
[0.9999995009985435, -0.000999000832833342]
[0.9999396214263468, -0.010988789821184267]
[0.9999396214263468, -0.010988789821184267]
[0.9968129466114029, -0.07977436592568171]
[0.9968129466114029, -0.07977436592568171]
[0.9712456180254766, -0.23807970964822986]
[0.9712456180254766, -0.23807970964822986]
[0.8864142948207253, -0.4628927349710369]
[0.8864143009584838, -0.46289273817622417]
[0.7057801368177318, -0.7084309119416488]
[0.7057801514552834, -0.7084309266341763]
⋮
[0.06392629002627652, -0.9979767881959541]
[0.06392487602860723, -0.9979547137591]
[-0.7349058827090867, -0.6782391803349066]
[-0.7348709624987291, -0.6782069527321518]
[-0.9963993606271645, 0.08509220057406423]
[-0.9963732720449674, 0.08508997261712575]
[-0.5728055158917739, 0.8197420053777407]
[-0.5727817136718961, 0.8197079420195611]
[0.23614725051915167, 0.971749724329546]
[0.23613980636823712, 0.9717190915291082]
[0.8656881437544278, 0.5006247531455262]
[0.86567035212763, 0.5006144643032474]
using DifferentialEquations, StaticArrays
function harmonic!(u,p,t)
[u[2],-u[1]]
end
u0 = @MArray [1.0, 0.0]
function constantenergymanifold(resid,u,p,t)
resid[1] = u[2]^2 + u[1]^2 - 1.0
resid[2] = 0.0
end
constantenergy_cb = ManifoldProjection(constantenergymanifold)
prob = ODEProblem(harmonic!, u0, (0., 100.), nothing)
sol = solve(prob,Tsit5())
ManifoldProjection is required for this problem to occur, so I think the problem is with NLsolve.jl's MArray handling.
It occurs even with the oop interface if ManifoldProjection is there:
using DifferentialEquations, StaticArrays
function harmonic!(u,p,t)
[u[2],-u[1]]
end
u0 = @MArray [1.0, 0.0]
function constantenergymanifold(resid,u,p,t)
resid[1] = u[2]^2 + u[1]^2 - 1.0
resid[2] = 0.0
end
constantenergy_cb = ManifoldProjection(constantenergymanifold)
prob = ODEProblem(harmonic!, u0, (0., 100.), nothing)
sol = solve(prob,Tsit5(),callback = constantenergy_cb)
highly pointing towards NLsolve.
In fact,
using DifferentialEquations, StaticArrays
function harmonic!(u,p,t)
[u[2],-u[1]]
end
u0 = @MArray [1.0, 0.0]
function constantenergymanifold(resid,u,p,t)
resid[1] = u[2]^2 + u[1]^2 - 1.0
resid[2] = 0.0
end
constantenergy_cb = ManifoldProjection(constantenergymanifold)
prob = ODEProblem(harmonic!, u0, (0., 100.), nothing)
sol = solve(prob,Tsit5(),callback = constantenergy_cb)
sol.t[end-1] == sol.t[end]
all(isnan.(sol[end]))
the fact that the time points are the same at the end means that non-NaN values go into NLsolve.jl, and NaNs come out. So it's clearly a problem with NLsolve and there's probably an only-NLsolve.jl MWE hiding in here.
using NLsolve
function f!(F, x)
F[1] = (x[1]+3)*(x[2]^3-7)+18
F[2] = sin(x[2]*exp(x[1])-1)
end
function j!(J, x)
J[1, 1] = x[2]^3-7
J[1, 2] = 3*x[2]^2*(x[1]+3)
u = exp(x[1])*cos(x[2]*exp(x[1])-1)
J[2, 1] = x[2]*u
J[2, 2] = u
end
nlsolve(f!, j!, @MArray [ 0.1; 1.2])
MArrays always throw NaNs in NLsolve.jl, so that's your the issue right there. It's upstream so I'll mark this upstream.