Alternative error control on Rosenbrock23 algebraic equations

[ChrisRackauckas]

Adds a separate error term for Rosenbrock23 when the mass matrix is not I in order to account for the fact that the algebraic conditions are not error controlled. They are added via doing error control on the absolute tolerance of the residual in the algebraic conditions. There isn't a sensible definition of relative tolerance here so leaving it at that.

Note that this technique is not applicable to the Rodas / stiffly stable methods since they always exactly satisfy the algebraic conditions by design, i.e. this term is just zero. Those methods will need to use a different strategy like residual control.

Handles part of #2054.

Fixes #2079

[ChrisRackauckas]

using ModelingToolkit, Plots

@variables t
D = Differential(t)

@mtkmodel Foo begin
    @parameters begin
        ω₁
        ω₂
        l
    end
    @variables begin
        x(t) # differential variable
        y(t) # differential variable
        z(t) # algebraic variable
    end
    @equations begin
        D(x) ~  ω₁ * 2pi * y    # solution: x(t) = sin(2πω₁t)
        D(y) ~ -ω₁ * 2pi * x    # solution: y(t) = cos(2πω₁t)
        z^l ~ sin(ω₂ * 2pi * t) # solution: z(t) = sin(2πω₂t) when l == 1.0
    end
end

@mtkbuild foo = Foo();

using OrdinaryDiffEq: solve, Rosenbrock23, Rodas4, FBDF

prob = ODEProblem(foo,
                  [foo.x => 1.0, foo.y => 0.0, foo.z => 0.0],
                  (0.0, 1.0),
                  [foo.ω₁ => 1.0, foo.ω₂ => 1000.0, foo.l => 1.0])

sol0 = solve(prob, Rosenbrock23(), maxiters = 1e7)
sol1 = solve(prob, Rodas4())
sol2 = solve(prob, FBDF())

plots = []
for (name, sol) in [
    ("Rosenbrock23", sol0),
    ("Rodas4", sol1),
    ("FBDF", sol2),
]
    z = sol(sol.t, idxs=foo.z)
    plt = plot(z, layout=(2,1), legend=:none, title="$(name) (t=0..1)")
    plot!(plt[2], z[1:findfirst(sol.t .>= 2e-3)], xlim=(0,2.1e-3), legend=:none, title="$(name) (t=0..2e-3)")
    push!(plots, plt)
end
plot(plots...)

Before:

After:

[oscardssmith]

is it reasonable that Rosenbrock23 now takes 1000x more steps than the other solvers?

[ChrisRackauckas]

Something is up with that test, I think I need to bring it local and check. It wouldn't be that unreasonable though since it's not B-stable.