Voltage divider gives wrong answer with Rosenbrock23

[oxinabox]

I haven't yet got a MWE without ModellingToolkit standard library. but I am pretty sure this isn't MSL, as this is a reproduction of the same failure on the same circuit when implemented in Cedar, which doesn't use MSL

This is a classic circuit, with the well known solution that the voltage across each of 2 equal resistors in series is half the voltage across both of them.

using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkitStandardLibrary.Blocks: Step,
    Constant, Sine, Cosine, ExpSine, Ramp,
    Square, Triangular
using ModelingToolkitStandardLibrary.Blocks: square, triangular
using OrdinaryDiffEq: ReturnCode.Success

@parameters t
@testset "voltage divider" begin
    @named source = Sine(offset = 0, amplitude = 1.0, frequency = 1e3, start_time = 0.5,phase = 0)
    @named voltage = Voltage()
    @named R1 = Resistor(R = 1e3)
    @named R2 = Resistor(R = 1e3)
    @named ground = Ground()

    connections = [connect(source.output, voltage.V)
        connect(voltage.p, R1.p)
        connect(R1.n, R2.p)
        connect(R2.n, voltage.n, ground.g)]

    @named model = ODESystem(connections, t,
        systems = [R1, R2, source, voltage, ground])
    sys = structural_simplify(model)
    prob = ODEProblem(sys, Pair[R2.i => 0.0], (0, 2.0))
    sol = solve(prob,  Rosenbrock23())
    @test sol.retcode == Success
    @test sol[voltage.p.v] ≈ 2*sol[R2.v]
end

output is

voltage divider: Test Failed at /home/oxinabox/JuliaEnvs/ModelingToolkitStandardLibrary.jl/test/Electrical/analog.jl:29
  Expression: sol[voltage.p.v] ≈ 2 * sol[R2.v]
   Evaluated: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6422526531740537, -5.09502586896226e-13] ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.9082824125929643, 3029.774147437094]

WE can see at the end the voltage goes completely wild 3029.7 is not a voltage that can occur given the voltage source is a amplitude 1 sin wave.

This does not error with Rodas4, and Tsit5 gives an appropriate error.

I suspect that this should also error, and this is some issue with Rosenbrock that it doesn't actually support purely algebraic equations.

[oxinabox]

Got a more minimal example

using OrdinaryDiffEq
function f!(out, u, _, t)
    out[1] = u[1]  - sin(1e8*t)
end

mprob = ODEProblem(ODEFunction(f!, mass_matrix=[0.0;;]), [0.0], (0, 2.0))
msol2 = solve(mprob,  Rosenbrock23())

We see something that i think is fairly nonsense from Rosenbrock. Since we know that u[1]=sin(1e8*t) is the solution so for every value of t we should see -1 <= u[1] <= 1. Not numbers like -1e7 Rodas4 gives perfectly good answers for this OTOH

julia> msol2 = solve(mprob,  Rosenbrock23())
retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 9-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.0011109999999999998
 0.011110999999999996
 0.11111099999999996
 1.1111109999999995
 2.0
u: 9-element Vector{Vector{Float64}}:
 [0.0]
 [-71.60648190259636]
 [-608.378199027157]
 [-6393.650939103113]
 [51590.080298312256]
 [-532821.4537692027]
 [-6.618206837399774e6]
 [7.051603353318885e7]
 [-1.8435495481137395e7]

[ChrisRackauckas]

What if you set a much smaller maxdt?

[oxinabox]

I still get different answers to Rodas4 and that still look unreasonable (but less so) at 1e-5. And if I go down to 1e-6 it errors

julia> msol2 = solve(mprob,  Rosenbrock23(), dtmax=1e-5)
retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 200002-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 2.1000000000000002e-5
 3.1e-5
 4.1e-5
 5.1e-5
 6.1e-5
 7.1e-5
 8.1e-5
 9.1e-5
 0.000101
 0.000111
 0.000121
 0.000131
 0.000141
 0.000151
 0.000161
 0.000171
 0.000181
 0.000191
 0.000201
 0.000211
 0.000221
 0.000231
 0.000241
 ⋮
 1.9997510000046332
 1.9997610000046333
 1.9997710000046334
 1.9997810000046334
 1.9997910000046335
 1.9998010000046336
 1.9998110000046336
 1.9998210000046337
 1.9998310000046338
 1.9998410000046338
 1.999851000004634
 1.999861000004634
 1.999871000004634
 1.999881000004634
 1.9998910000046342
 1.9999010000046342
 1.9999110000046343
 1.9999210000046344
 1.9999310000046344
 1.9999410000046345
 1.9999510000046345
 1.9999610000046346
 1.9999710000046347
 1.9999810000046347
 1.9999910000046348
 2.0
u: 200002-element Vector{Vector{Float64}}:
 [0.0]
 [-71.60648190259636]
 [-608.378199027157]
 [-642.7489774653449]
 [-114.55895344035164]
 [513.8978606321156]
 [692.569761780467]
 [265.0756251618269]
 [-394.423791329092]
 [-708.7070000312888]
 [-402.70018474761207]
 [255.76668419037247]
 [690.375847949447]
 [520.7391791360687]
 [-104.67021084760411]
 [-638.4678521198861]
 [-613.4517109841029]
 [-51.51696102449541]
 [555.5075890755888]
 [676.3286506584927]
 [205.19857457772727]
 [-445.5298810041105]
 [-706.311940455724]
 [-348.9002322890573]
 [313.88355975105]
 [701.9433250803237]
 ⋮
 [708.4770946617418]
 [382.9440902082443]
 [-277.7576066589109]
 [-695.3542227393542]
 [-504.347724299085]
 [128.08500816175558]
 [648.4123899072378]
 [601.2220929110389]
 [27.817060606508694]
 [-569.9346272130174]
 [-668.855672018441]
 [-182.36623023532903]
 [463.7377678057165]
 [703.9590652347118]
 [328.0459298992676]
 [-334.9867497711738]
 [-704.8249889896883]
 [-457.7709027380784]
 [189.94343412048536]
 [671.4113355574607]
 [565.2319189343934]
 [-35.662126636769315]
 [-605.3431866074901]
 [-645.2025576711708]
 [-120.3536293154349]
 [459.4856596349784]

these values are also obviously wrong since they are outside of [-1,1]

[topolarity]

1.6 MHz frequency over 2 seconds is a little extreme and requires bumping maxiters, so I knocked it down to 1 kHz.

Rosenbrock23 takes exponentially increasing timesteps and fails to converge to a good u, as reported:

# Rosenbrock23
     sol.t: [0.0, 1.0e-6,     1.1e-5,   0.000111,  0.001111,  0.011111, 0.111111,   1.111111,   2.0]
  sol.u[1]: [0.0, 0.00628317, 0.069085, 0.66682,  -7.14876,  -33.4121, -339.90162, -3404.7964, -3029.7741]
 f(t,u)[1]: [0.0, 2.3698e-8,  2.55e-5,  0.0246,   -7.791,    -34.054,  -340.544,   -3405.439,  -3029.77]

Rodas4 does successfully solve the nonlinear equation at each timestep but still takes very unreasonable timesteps (in fact, exactly the same timesteps as Rosenbrock23):

# Rodas4
     sol.t: [0.0, 1.0e-6,   1.1e-5,  0.000111, 0.001111, 0.011111, 0.111111, 1.111111, 2.0]
  sol.u[1]: [0.0, 0.006283, 0.06906,  0.64225, 0.64225,  0.6423,   0.6423,   0.6423,  -1.286e-12]
 f(t,u)[1]: [0.0, 0.0,      0.0,     0.0,      0.0,      0.0,      0.0,      0.0,      0.0]

FBDF, in contrast, takes 118668 timesteps (~118 samples per period) and produces a good result (maximum(|f(t,u)[1]|) = 2.6937e-6):