Open ArnoStrouwen opened 1 month ago
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.
But it's still strange that the solve
that uses the Jacobian fails.
But is y
not a function of z
and x
here?
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)
Actually, reviewing this.
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
https://docs.sciml.ai/ModelingToolkit/stable/tutorials/nonlinear/
At first glance, the generated jac looks overly simplistic?