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