ChrisRackauckas
commented
6 months ago using OrdinaryDiffEq
dt = 1e-6
abstol = 1e-7
d=1
k=1000
F = 100
autodiff=false
function f2(du,u,p,t)
F, k, d = p
x, dx, ddx = u
du[1] = dx
du[2] = ddx
du[3] = (d*dx + k*x^2) - (F)
end
fmm2 = ODEFunction(f2; mass_matrix=[1 0 0;0 1 0;0 0 0])
prob2 = ODEProblem(fmm2, [0.0, F/d, 0.0], (0.0, 0.01), [F, k, d])
ref1 = solve(prob2, Rodas5P(;autodiff); abstol, dt, initializealg=NoInit())
plot(ref1)
ref2 = solve(prob2, ImplicitEuler(;autodiff, nlsolve=NLNewton(check_div = false, always_new=true)); abstol, dt, adaptive=false, initializealg=NoInit())
plot(ref2)The solution is fine with that. If adaptivity is turned off then the nonlinear solver needs to be tweaked in order to not easily give up. It easily gives up on purpose because the quickest thing to do is just change dt. It's a common strategy that all stiff ODE solvers do, and it's also why most software doesn't allow you to set adaptive=false.
We should when adaptive=false set those two options. But instead of doing that, I want to just complete https://github.com/SciML/OrdinaryDiffEq.jl/issues/1570 and https://github.com/SciML/NonlinearSolve.jl/issues/241, in which case I would propose a very different handling. Since that's scheduled for the next two weeks, no reason to do anything here and just finish the real answer.
bradcarman
mtiller-jh
Bug and MWE
The following problem does not solve (i.e. gives
Warning: Newton steps could not converge and algorithm is not adaptive. Use a lower dt.)Expected behavior
This problem solves without issue using a clean implementation of the Implicit/Backwards Euler method. See https://github.com/bradcarman/SimpleEuler.jl/blob/main/test/runtests.jl for more information.
Note: this problem also solves correctly using
f0andf1versions of the problem, as also seen in the SimpleEuler.jl runtests.jl file. Adding additional derivatives to be calculated should not be causing theImplicitEulermethod to fail.Additionally the Modelica problem solves without issue and gives the same result: link
Here is a plot of
ddxfrom ref0, ref1, sol0, sol1, sol2, and sol3. * -ddxcalculated manuallyHere is the modelica result:
Environment