Open hersle opened 2 weeks ago
The first problem is that I want to get rid of the evalsol()
workaround that hardcodes idxs=y[0]
. Ideally, I'd like to use callable parameters, but to start with I try
@register_symbolic evalsol(sol::ODESolution, t::AbstractFloat, idx::Num)
evalsol(sol::ODESolution, t::AbstractFloat, idx::Num) = sol(t, idxs=idx)
# 1st order system: replace 0th order ODE by its solution
@parameters sol0::ODESolution
@named sys1 = ODESystem([y[0] ~ evalsol(sol0, t, y[0]); eqs1], t; defaults = ics1)
ssys1 = structural_simplify(sys1)
Now the problem is that the 0th order equation in ssys1
is no longer an observed equation, which it should be:
equations(ssys1)
2-element Vector{Equation}:
Differential(t)((y(t))[1]) ~ 2(y(t))[0]
0 ~ -(y(t))[0] + evalsol(sol0, t, (y(t))[0]) # <-- y[0] should be observed!
That makes ssys1
an unbalanced/unsolvable system (in the same way as above). Is there a way around this?
Here is a simple ODE system (from some perturbation theory toy problem; the real problem is bigger):
I can easily solve the whole system "simultaneously":
But now I want to exploit this property of the system:
y[0]
is independent ofy[1]
(buty[1]
depends ony[0]
). This always happens in perturbation theory. I want to solve for onlyy[0]
first, and then use thatODESolution
to solve for onlyy[1]
. I split the "simultaneous" system into two "sequential" ones (my dream is to make a function that does this for a general system, i.e. a block lower triangular transformation):I now solve them "sequentially", and everything works!
But there are some problems that block me from doing this generally (next post).