Broken Modeling Nonlinear Systems tutorial when using analytical jac

[ArnoStrouwen]

https://docs.sciml.ai/ModelingToolkit/stable/tutorials/nonlinear/

At first glance, the generated jac looks overly simplistic?

julia> prob.f.jac.jac_iip
RuntimeGeneratedFunction(#=in ModelingToolkit=#, #=using ModelingToolkit=#, :((ˍ₋out, ˍ₋arg1, ˍ₋arg2)->begin
          #= C:\Users\arno\.julia\packages\SymbolicUtils\dtCid\src\code.jl:373 =#
          #= C:\Users\arno\.julia\packages\SymbolicUtils\dtCid\src\code.jl:374 =#
          #= C:\Users\arno\.julia\packages\SymbolicUtils\dtCid\src\code.jl:375 =#
          begin
              begin
                  y = (*)(ˍ₋arg1[1], (+)((*)(-1, ˍ₋arg1[2]), ˍ₋arg2[1]))
                  begin
                      #= C:\Users\arno\.julia\packages\Symbolics\Nk48t\src\build_function.jl:546 =#
                      #= C:\Users\arno\.julia\packages\SymbolicUtils\dtCid\src\code.jl:422 =# @inbounds begin
                              #= C:\Users\arno\.julia\packages\SymbolicUtils\dtCid\src\code.jl:418 =#
                              ˍ₋out[1] = (*)(-1, ˍ₋arg2[3])
                              ˍ₋out[2] = y
                              ˍ₋out[3] = 0
                              ˍ₋out[4] = (*)(-1, ˍ₋arg2[2])
                              #= C:\Users\arno\.julia\packages\SymbolicUtils\dtCid\src\code.jl:420 =#
                              nothing
                          end
                  end
              end
          end
      end))

[hersle]

I don't think it's too simplistic. The simplified equations are

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

with Jacobian

julia> calculate_jacobian(ns)
2×2 Matrix{Num}:
 -σ   0
  y  -β

which seems to match the generated expressions.

[hersle]

But it's still strange that the solve that uses the Jacobian fails.

[ArnoStrouwen]

But is y not a function of z and x here?

[hersle]

Ah, of course, you're right. I think the correct Jacobian should be

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

@parameters σ ρ β
@variables x z
Dx(expr) = expand_derivatives(Differential(x)(expr))
Dz(expr) = expand_derivatives(Differential(z)(expr))

y = x*(-z + ρ)
f1 = (-x + y)*σ
f2 = x*y - z*β

J = [
    Dx(f1) Dz(f1)
    Dx(f2) Dz(f2)
]

which gives

2×2 Matrix{Num}:
 (-1 - z + ρ)*σ        -x*σ
    2x*(-z + ρ)  -β - (x^2)

[ChrisRackauckas]

Fixed by #2863

[ChrisRackauckas]

Actually, reviewing this.

[ChrisRackauckas]

I thought I might as well check how ODEs is accomplishing this, and it's a simpler solution through full_equations which expands the observed. https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/diffeqs/abstractodesystem.jl#L28. We should probably change the solution on NonlinearSystem for that.

We should also check the OptimizationSystem gradient for the same effect. https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/optimization/optimizationsystem.jl#L130