No, second order is not really all that high. You may still experience numerical difficulties due to poor scaling though, even for this simple system. Consider what happens if you increase w1, for w1=1e6, the solver exits with
┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=1.9999999999999998, and step error estimate = 8.930759781475992. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/s9wrq/src/integrator_interface.jl:599
retcode: DtLessThanMin
Already for w1=1000, you can see part of the problem in the solution:
model, u, y = model_bessel(2, 100.0, w1=1e3)
prob = ODEProblem(model, [], tspan)
sol = solve(prob, Tsit5())
plot(sol)
The derivative state becomes enormous, much larger that the position state.
Compare this to a balanced representation:
using ControlSystemsMTK, ControlSystemsBase
function balanced_model_bessel(step_time, step_size; w1=3.14, y0=0.0, yd0=0.0)
println("Defining the model...")
@variables t u(t)=0 y(t)=y0 yd(t)=yd0
@parameters w=w1
H = tf(w1^2, [0.618, 1.3617w1, w1^2])
Hss = ss(H) # This performs balancing
@named filter = ODESystem(Hss)
eqs = [filter.input.u ~ u
u ~ ifelse(t > step_time, step_size, 0.0)]
@named sys = ODESystem(eqs, t, systems=[filter])
println("Model defined...")
sys = structural_simplify(sys)
println("Model simplified!")
sys, u
end
model, u = balanced_model_bessel(2, 100.0, w1=1e3)
prob = ODEProblem(model, [], tspan)
sol = solve(prob, Tsit5())
@show length(sol.t)
plot(sol)
Here, the two state variables are of comparable magnitude, making the job for the solver much easier.
To plot the correctly scaled filter output in this case, you’d have to plot plot(sol, idxs=model.filter₊output₊u), since the output is scaled w.r.t. the state.
