Change allow_symbolic to default to true

[ChrisRackauckas]

From the longer conversation. The reason for not defaulting to it before was a scare that the eliminated expression may be zero. This was found before, leading to NaNs in the resulting system evaluations. However, such zeros would also be problematic to the solver as well, since leaving the uneliminated term in leads to the structural jacobian not matching the true jacobian, which can make the structural DAE index different from the true index. Because of this phenomena, it is no less safe to eliminate the extra variables. But it can lead to some numerical improvements. For example, p * y ~ p * z is trivial at p = 0, but y ~ z is not, and therefore eliminating the p is more numerically robust.

[hersle]

This is also a nicer default for structural simplification after independent variable change, provided that d(old_iv)/d(new_iv) is nonzero, which it has to be for the transformation to be well-defined.

[AayushSabharwal]

I looked at the failing tests. allow_symbolic = true might be problematic for nonlinear systems.

This system in the tests (lorenz):

    @variables x y z
    @parameters σ ρ β
    eqs = [0 ~ σ * (y - x)
           0 ~ x * (ρ - z) - y
           0 ~ x * y - β * z]
    guesses = [x => 1.0, z => 0.0]
    ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
    @mtkbuild ns = NonlinearSystem(eqs)

Ends up with

julia> equations(ns)
1-element Vector{Equation}:
 0 ~ x*y - z*β

julia> observed(ns)
2-element Vector{Equation}:
 y ~ x
 z ~ (-y + x*ρ) / x

julia> full_equations(ns)
1-element Vector{Equation}:
 0 ~ ((x - x*ρ)*β) / x + x^2

Whereas with allow_symbolic = false

julia> equations(ns2)
2-element Vector{Equation}:
 0 ~ (-x + y)*σ
 0 ~ x*y - z*β

julia> full_equations(ns2)
2-element Vector{Equation}:
 0 ~ (-x + x*(-z + ρ))*σ
 0 ~ -z*β + (x^2)*(-z + ρ)

In the first case, the jacobian at x => 0, z => 0 is NaN whereas in the second case it isn't

julia> prob = NonlinearProblem(ns, [x => 0.0], ps; jac=true)
NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 1-element Vector{Float64}:
 0.0

julia> prob2 = NonlinearProblem(ns2, [x => 0.0, z => 0.0], ps; jac=true)
NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
 0.0
 0.0

julia> prob.f.jac(prob.u0, prob.p)
1×1 Matrix{Float64}:
 NaN

julia> prob2.f.jac(prob2.u0, prob2.p)
2×2 Matrix{Float64}:
 250.0  -0.0
   0.0  -2.66667

[ChrisRackauckas]

The point though is that those cases where you get a NaN, you would have had numerical issues anyways. So the correct behavior would be to take the symbolic expressions found through this process and warn/error before solving if they are detected as being zero. That would be more clear to users anyways.

But we should make this a bit smarter:

julia> full_equations(ns)
1-element Vector{Equation}:
 0 ~ ((x - x*ρ)*β) / x + x^2

we should just multiply both sides by the denominator?

julia> full_equations(ns)
1-element Vector{Equation}:
 0 ~ ((x - x*ρ)*β)

is the same, and the best answer.

[baggepinnen]

Is the denominator there really x + x^2? If so, the printing is very unfortunate and should instead be / (x + x^2)

[ChrisRackauckas]

yes that's bad printing

[AayushSabharwal]

No, the denominator is just x? x^2 is different term

[hersle]

Maybe it should simplify like this? However, this assumes x ≠ 0.

 0 ~ ((x - x*ρ)*β) / x + x^2    --->     0 ~ (1 - ρ)*β + x^2

[AayushSabharwal]

Do we add the denominators as assertions and remove them? And then add the infrastructure to validate assertions on initialization.

[ChrisRackauckas]

What's left here?

[AayushSabharwal]

Alright so

• reference tests need an update
• need to propagate this to SDEs
• the signaling thing? Not sure how we should expose that

[AayushSabharwal]

Also Catalyst and MTKStdlib

[AayushSabharwal]

https://github.com/SciML/ModelingToolkit.jl/actions/runs/14399916432/job/40383199248?pr=3481#step:6:1405 this failing test in the allow_algebraic PR corresponds to https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/main/test/multi_domain.jl#L179, which is a purely algebraic system. The initial conditions are such that one of the observed variables (sys.source.i) has zero in the denominator at t=0. I’m not fully sure how we want to proceed with this. The actual value of sys.source.i at t = 0 is in fact zero.

Catalyst also seems to have a case where allow_algebraic reduces the number of equations in a nonlinear system from 6 to 1, but the resulting system now stalls with the default algorithm and any other algorithm I try.