Silent and wrong derivatives of observed variables

[hersle]

Derivatives of observed variables calculated with ODESolution(t, Val{1}, ...) are incorrect.
On v9.12.2 (after #2574), the simple test

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

@variables x(t) y(t)
@named sys = ODESystem([D(x) ~ 1, y ~ x^2], t)
sys = structural_simplify(sys) # move y to observed variables
prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0)) # analytical solution: x = t, y = t^2
sol = solve(prob)

y′_obs(t) = sol(t, Val{1}, idxs=y) # request derivative of observed variable y
y′_anal(t) = 2 * t # analytical solution y′ = 2*x*x′ = 2*t
@test y′_obs(1.0) ≈ y′_anal(1.0)

fails with

Test Failed at bug.jl:14
  Expression: y′_obs(1.0) ≈ y′_anal(1.0)
   Evaluated: 0.9999999999999751 ≈ 2.0

Is this the general case mentioned in #2574?
Having this work generally would be both very convenient and a natural extension to the solution calling interface!

Also, the current behavior is to return the wrong derivative with no error or warning.
I think this is severe and likely to cause hard-to-find bugs in user code.
Is it possible to identify the cases when y is an observed variable and the methods of #2574 do not apply?
In those cases, for example, could ODESolution(t, Val{N >= 1}, idxs=y) issue an error or warning?

[ChrisRackauckas]

Interesting, not entirely sure what's going on here. We can take a look, it's fixable.

[AayushSabharwal]

This is because interpolation produces the timeseries for D(x) (derivative of the only state variable), and then indexes it for y. Since y is observed, it generates a function that returns x^2.

[AayushSabharwal]

I feel the best we can do here right now is error if any of the indexes is an observed quantity and deriv > 0. SciMLBase doesn't know about derivatives and so can't ask for a function for D(y), assuming MTK was made to be able to generate this function (it can't currently, but that's easy). Eventually if we get interpolation into SII we can expose an API that solves this issue

[ChrisRackauckas]

That's fine, get an error setup for now and then we can fix it later.

[hersle]

Sounds good 👍 I think it is better than returning a wrong result
Maybe the error can include a helpful tip, like splining y yourself with SciML/DataInterpolations.jl and calling DataInterpolations.derivative on it, or something similar? That does at least cover my use case.

As always, thank you for your responsiveness 😃

[hersle]

Do you have some quick pointers on how to compute y′ from the ODESolution in the "symbolic-manual" way? Suppose I do

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

@variables x(t) y(t)
@named sys = ODESystem([D(x) ~ 1, y ~ x^2], t)
sys = structural_simplify(sys) # move y to observed variables
prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0)) # analytical solution: x = t, y = t^2
sol = solve(prob)

y_sym = Dict(eq.lhs => eq.rhs for eq in observed(sys))[y] # find expression for y in terms of other quantities
y′_sym = expand_derivatives(D(y_sym)) # expression for y′ in terms of other quantities
y′_obs(t) = NaN # TODO: how to compute this from ODESolution sol?
y′_anal(t) = 2 * t # analytical solution y′ = 2*x*x′ = 2*t

@test y′_obs(1.0) ≈ y′_anal(1.0)

Now y′_sym is 2x(t)*Differential(t)(x(t)). Is there an easy user-facing interface to construct the function y′_obs(t) that evaluates this expression numerically using sol?