Change independent variable of ODE systems #3437
Conversation
|
Really cool stuff!
Yes all of the time. I'll be happy to see this completed. |
|
I think this is ready for a look now. |
There was a problem hiding this comment.
This is great! I have a version of this that I use and was hoping it would make its way to an MTK feature. I see that flattening is a requirement, though. I was able to get around this by essentially applying what you have recursively through each subsystem of the hierarchical system. Would that be possible in a follow-up?
My implementation did not feel robust enough to open a PR with (did not have the chain rule in there). Would be happy to try and contribute here though.
Thanks!
|
|
||
| # 2) Find e.g. du/dt in equations, then calculate e.g. d²u/dt², ... | ||
| eqs = [get_eqs(sys); eqs] # all equations (system-defined + user-provided) we may use | ||
| idxs = findall(eq -> isequal(eq.lhs, D1(iv2_of_iv1)), eqs) |
There was a problem hiding this comment.
will this catch instances of c * D(x) ~ something_rhs?
There was a problem hiding this comment.
Not currently. I also thought about this and would love to fix it! Very open to suggestions on a how to make a more general solution for this.
There was a problem hiding this comment.
I essentially have a function get_de_arg(eq::Equation) which assumes that only the lhs OR rhs has a differential, picks the side with the differential using has_differential. has_differential uses mtk/symbolics isdifferential anf if true returns early. If not, is checks hasproperty(f, :arguments), which if yes, cycles through them and again checks has_differential returning true if so. else the top-level function of the recursive call returns false...
I think I modified something I found inside MTK/symbolics to do this?
But anyways, that is how I "get" the de arg for an equation to check against my desired iv.
There was a problem hiding this comment.
Thanks! That would be an improvement. Would that be equivalent to using Symbolics.solve_something for the derivative on the equation? Is there also a way to do the independent variable transformation in the general case, where there can be a nonlinear relation involving the derivative? Maybe @ChrisRackauckas
There was a problem hiding this comment.
will defer to @ChrisRackauckas ! Solution I described I implemented months ago when I was learning MTK so it is definitely not the most general...
There was a problem hiding this comment.
Not currently. I also thought about this and would love to fix it!
I hope there is an error thrown if anything that isn't handled correctly is encountered so that you do not get an incorrect transformation back? :)
There was a problem hiding this comment.
I probably won't be able to get to this until the weekend, but can try to give this a comprehensive review and see what I can upstream from my private implementation to try and cover the corner cases I've hit.
There was a problem hiding this comment.
@baggepinnen The next line in the code errors if a transformation of the wanted form is not specified.
Definitely. I want this to work too, but have just focused on getting the transformation itself right first. I see one non-trivial aspect of hierarchical systems. If the new independent variable is a dependent variable in a subsystem of the original system, is it OK that it will effectively be made global when it is turned into an independent variable in the new system and all its subsystems? Also, if the user wants the old independent variable as a new dependent variable in the new system, should it always be put in the top-level system? Maybe I'm overthinking this, though. |
Nope! all valid and things that I don't really have great answers to...In my use case, my desired iv is indeed in a subsystem so I also have to recursively step through my subsystems to a) get to the subsystem holding my new iv's differential equation and b) track the namespaces as I need to re-namespace everything in the rhs of that equation as I recursively modify/chain rule all the diffeqs in the other subsystems with this equation. One more complicating piece of this was that my model(s) are fully defined with units so they needed to be tracked as well... I do place my original iv--now a "variable" with a diffeq--in the top-level system. TBH I am not quite sure of the consequences of these things...and I feel that I probably break some interface with symbolics maybe (to address things like the comment I left), but all-in-all my system with tens(ish) subsystems (of varying depth), couple hundred equations and connections, and some other complexities still simplifies to what I would expect. |
|
I have tried to add support for hierarchical systems. See the new example/tests/documentation. If you want to try it on a larger example, that would be very helpful!
This is another good point. This information is currently lost because I used |
|
|
||
| # 1) Utility that performs the chain rule on an expression, e.g. (d/dt)(f(t)) -> (d/dt)(f(u(t))) -> df(u(t))/du(t) * du(t)/dt | ||
| function chain_rule(ex) | ||
| vars = get_variables(ex) |
There was a problem hiding this comment.
It'd be better to use ModelingToolkit.vars instead of get_variables. iirc there are cases the latter doesn't handle.
There was a problem hiding this comment.
I now use vars(ex; op = Nothing) instead. It seems to work the same for the current tests and examples. I need op = Nothing so it doesn't count e.g. D(f(t)) as a variable; I just want to replace f(t) in that case.
| y(x) | ||
| ``` | ||
| """ | ||
| function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, simplify = true, fold = false, kwargs...) |
There was a problem hiding this comment.
Does this support DDEs? If it does, should test it. If it doesn't, the function should check is_dde(sys) and error.
There was a problem hiding this comment.
No. I think that would require more effort, so I will forbid them to limit the scope for now. Added a test to see that it errors.
| eqs = [eqs; get_eqs(sys)] # all equations (system-defined + user-provided) we may use | ||
| idxs = findall(eq -> isequal(eq.lhs, D1(iv2_of_iv1)), eqs) | ||
| if length(idxs) != 1 | ||
| error("Exactly one equation for $D1($iv2_of_iv1) was not specified!") |
There was a problem hiding this comment.
It might be helpful to print eqs[idxs] as part of the error message, so the user knows which equations we're looking at.
There was a problem hiding this comment.
Good idea, I updated it to do that.
| # 3) If requested, insert extra dummy equations for e.g. du/dt, d²u/dt², ... | ||
| # Otherwise, replace all these derivatives by their explicit expressions | ||
| if dummies | ||
| div2name = Symbol(iv2name, :_, iv1name) # e.g. :u_t # TODO: customize |
There was a problem hiding this comment.
An edge case, but what if a variable with the name u_t already exists in the system? A better approach would be to use default_toterm. Basically something like div2name = nameof(default_toterm(D1(iv2_of_iv1))) and a double derivative for ddiv2name.
There was a problem hiding this comment.
Is that what gives the "hidden" unicode underscore? I am a little skeptical, because I would like the "dummy" variables to be exposed to the user (and not only something internal).
I see what you mean, though. And I suppose it is confusing if change_independent_variable and structural_simplify create variables with different but indistinguishable underscores.
There was a problem hiding this comment.
I went with your suggestion for now. It is good that it is consistent with how other dummy derivatives in MTK works, and that it prevents name collisions unless the user has a secret love relationship with U+02CD 😅 But I think it's bad that it is obscure for the user to access.
Do you see any other possible solutions to this?
There was a problem hiding this comment.
The new variables can still be referred to as Differential(t)(u(t)), so they aren't inaccessible. If u is the new independent variable structural_simplify will create dummy derivative variables suffixed with _u, so it won't conflict.
| iv1 = get_iv(sys) # e.g. t | ||
|
|
||
| if !iscomplete(sys) | ||
| error("System $(nameof(sys)) is incomplete. Complete it first!") |
There was a problem hiding this comment.
Why require complete? If this handles unflattened systems (as the recursive nature of transform would suggest) then complete shouldn't be necessary.
There was a problem hiding this comment.
You are right, transform can handle it. I updated the examples to not use complete. When run on incomplete systems, it can be a little subtle to keep track of the differences between the old/new variables before/after structural simplification, so I added a warning note in the documentation about that.
|
This also needs a format. And could you add a test using a callable parameter? Including cases where it is called with and without dependent variables. |
|
Thanks for the comprehensive review and useful comments! I will get around to adding a test with a callable parameter and/or registered function soon. I think those may be a little subtle because I think they should be excepted from the chain rule machinery. I will do the format when everything is ready. |
…d nonlinear equations
Changing the independent variable of an ODE system is very useful in some applications:
This is often done with lots of error-prone manual chain rule calculations etc. I want to fully automate this process.
Checklist
contributor guidelines, in particular the SciML Style Guide and
COLPRAC.
TODO before merge