ImplicitEuler fails to solve a nonlinear spring-damper problem #1010
Comments
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 We should when |
|
Thanks Chris, that sounds good. A couple points:
I would assume that if using OrdinaryDiffEq
using SimpleEuler: BackwardEuler
using Plots
dt = 1e-4
abstol = 1e-3
d=100
k=1000
F = 100
autodiff=true
always_new=false
function f3(du,u,p,t)
F, k, d = p
x, dx, ddx, dddx = u
du[1] = dx
du[2] = ddx
du[3] = dddx
du[4] = (d*dx + k*x^2) - (F)
end
fmm3 = ODEFunction(f3; mass_matrix=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 0])
prob3 = ODEProblem(fmm3, [0.0, F/d, 0.0, 0.0], (0.0, 2), [F, k, d])
ref3 = solve(prob3, ImplicitEuler(;autodiff, nlsolve=NLNewton(;always_new, check_div=false)); abstol) #always fails: DtLessThanMin
ref3 = solve(prob3, ImplicitEuler(;autodiff, nlsolve=NLNewton(;always_new, check_div=false)); abstol, dt, adaptive=false) #works
sol3 = solve(prob3, BackwardEuler(;autodiff, always_new); abstol, dt, adaptive=false) #works
plot(ref3; idxs=3)
plot!(sol3; idxs=3)
|
Yes ImplicitEuler is first order so solving to a low tolerance is pretty infeasible and requires a really small dt.
That's kind of a misnomer since it "works" because it's ignoring the error and not actually solving to tolerance. Only the nonlinear solver iterations are bounded, not the truncation error. |
|
Just as a point of reference, I'd just like to present a different example as a different perspective on this issue. The problem that Brad is proposing is effectively: From this equation (and this equation alone), we can determine what But what I'd just like to introduce mostly as a thought experiment is the idea that imagine there was a closed form solution here (that we could derive from the equation above) such that: If that were the case and we could skip the numerical solution altogether, then what are the implications of adding the following relations: Because we know The key thing I want to point out here is that we aren't actually integrating anything, we are only differentiating. Now, imagine you plugged the equations above into The fact that @bradcarman can solve this by doing a backward Euler substitution seems to me just evidence that such a substitution is actually being used to differentiate in his case, not to integrate. In other words, it leads to this: ...which is equivalent to: My point here is that substituting backward Euler into an equation system isn't necessarily solving any differential equations. I'm not saying it doesn't "work", I'm just pointing out that it isn't necessarily accomplishing the same thing that |
|
As of DifferentialEquations v7.12.0 this is now working... Works with both Code with analytical comparison... using DifferentialEquations
# F = d*dx + k*x^2
function f(du,u,p,t)
F, k, d = p
x, dx, ddx = u
du[1] = dx #D(x) ~ dx
du[2] = ddx #D(dx) ~ ddx
du[3] = (F) - (d*dx + k*x^2)
end
abstol = 1e-3
d=1
k=1000
F = 100
odef = ODEFunction(f; mass_matrix=[1 0 0;0 1 0;0 0 0])
prob = ODEProblem(odef, [0.0, F/d, 0], (0.0, 0.01), [F, k, d])
sol = solve(prob; abstol) # Success
sol′ = solve(prob, ImplicitEuler(); abstol) # Success
c₁ = 0.0
t = 0:0.0001:0.01
x(t) = (sqrt(F)/sqrt(k))*tanh((c₁*d*sqrt(F)*sqrt(k)+ sqrt(F)*sqrt(k)*t)/d)
dx(t) = ForwardDiff.derivative(x, t)
ddx(t) = ForwardDiff.derivative(dx, t)
fig = Figure(resolution=(1600,400))
ax = Axis(fig[1,1]; ylabel="x(t)", xlabel="t")
lines!(ax, t, x.(t))
scatter!(ax, sol.t, sol[1,:])
ax = Axis(fig[1,2]; ylabel="dx(t)", xlabel="t")
lines!(ax, t, dx.(t))
scatter!(ax, sol.t, sol[2,:])
ax = Axis(fig[1,3]; ylabel="ddx(t)", xlabel="t")
lines!(ax, t, ddx.(t))
scatter!(ax, sol.t, sol[3,:])
fig
|
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
The text was updated successfully, but these errors were encountered: