Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Capture observed dependence on unknowns in simplified nonlinear systems #2863

Merged
merged 12 commits into from
Jul 16, 2024

Conversation

hersle
Copy link
Contributor

@hersle hersle commented Jul 16, 2024

This fixes the incorrect Jacobian in #2858.

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

@hersle
Copy link
Contributor Author

hersle commented Jul 16, 2024

Would it be better if, for example,

using ModelingToolkit
@variables x y z
@parameters σ ρ β
@mtkbuild ns = NonlinearSystem([
    0 ~ σ * (y - x)
    0 ~ x *- z) - y
    0 ~ x * y - β * z
])
J = calculate_jacobian(ns)

did not expand observed derivatives like

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

but rather kept them in a form like

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

?

# TODO: must do the same fix in e.g. calculate_hessian?
# TODO: rather keep observed derivatives unexpanded, like "Differential(obs)(expr)"?
obs = map(eq -> eq.lhs => eq.rhs, observed(sys))
rhs = map(eq -> substitute(eq.rhs, obs), equations(sys))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This won't catch a chain of observed, it would need a fixed point sub I think

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like in the most recent commit? Can you think of a simple example on to test?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

z ~ y, y ~ x, x^2 - 2z` should hit it, since if you just plug in z you get y.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not quite after structural simplification, it seems.
But I added 0 ~ x^2 + 2*z + y, z ~ y, y ~ x, which seems to do the trick.

@ChrisRackauckas ChrisRackauckas merged commit 0cc516e into SciML:master Jul 16, 2024
21 of 22 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants