Generalize symbolic indexing of ODE/DDE to multivariate variables

[hersle]

This is one possible fix to #3737.
It generalizes a line that assumed DDE variables are like $x(t-\tau)$ to permit multivariate $x(t-\tau, a_2, a_3, \ldots)$.
Indexing ODEs with multivariate variables works "for free" when the DDE line no longer errors.

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]

Any input on this? @ChrisRackauckas / @AayushSabharwal

[AayushSabharwal]

I'm a little on the fence about this. Is it just for UX?

[ChrisRackauckas]

That's not an ODE/DDE, so it should fail. It's not solvable with these techniques so it should fail at construction. Maybe we should throw a better error but an ODE solver cannot be directly applied to a PDE without proper handling of the other independent variable.

[hersle]

I totally see where you are coming from, but I think it is slightly more nuanced.

an ODE solver cannot be directly applied to a PDE

It can if the PDE is an ODE. In my application this is precisely a PDE that reduces to one ODE in $t$ for every $k$.

Is it just for UX?

The dependence on $(t, k)$ arises naturally after Fourier transformation away $\boldsymbol{x}$ in the PDE's original formulation in $(t, \boldsymbol{x})$.
It is critical for my application to distinguish ODE variables that depend on $(t, k)$ from those that only depend on $t$.
It is therefore conventional in my field to write $(t, k)$ where $t$ is the independent variable and $k$ is a parameter.

It is 100% intentional and with the understanding that my ODEs are reduced PDEs that I write $(t, k)$ and use ODE solvers.
Permitting this notation would be a big boost for my package. Do you see a way we can resolve this? Thank you.

[ChrisRackauckas]

It can if the PDE is an ODE. In my application this is precisely a PDE that reduces to one ODE in t for every k.

That would require infinitely many ODE solves then.

The dependence on $(t, k)$ arises naturally after Fourier transformation away $\boldsymbol{x}$ in the PDE's original formulation in $(t, \boldsymbol{x})$. It is critical for my application to distinguish ODE variables that depend on $(t, k)$ from those that only depend on $t$. It is therefore conventional in my field to write $(t, k)$ where $t$ is the independent variable and $k$ is a parameter.

That is very different. For that you have an array variable for a system of equations. k is an integer, and the Fourier transform creates N equations. What you wrote is continuous k over the space of real numbers.

[hersle]

That would require infinitely many ODE solves then.

It sure does 😅 The ODEs are big and stiff, too. I use EnsembleProblem to parallellize them.

What you wrote is continuous k over the space of real numbers.

The Fourier transform is in the infinite-space limit, so $k$ is in fact continuous.
My strategy is to solve on a grid of $k$, and then interpolate in-between to arbitrary $k$ as the solution varies smoothly over $k$.
This is faster and more stable than solving for all $k$ in the same system.

[ChrisRackauckas]

My strategy is to solve on a grid of $k$, and then interpolate in-between to arbitrary $k$ as the solution varies smoothly over $k$. This is faster and more stable than solving for all $k$ in the same system.

This is how you're discretizing the PDE. The ODE solver doesn't have ways to discretize a PDE

[hersle]

This is how you're discretizing the PDE.

Yes, sure.
Notably, my equations are linear, so the ODE for each $k_i$ is independent of other $k_j \neq k_i$.
Thus every $k$ boils down to an ODE in $\tau$ where $k$ is a parameter, and each ODE is well-defined on its own.
I have a custom type that wraps multiple ODE problems/solutions to take care of the interpolation in $k$ etc.

[ChrisRackauckas]

Notably, my equations are linear, so the ODE for each $k_i$ is independent of other $k_j \neq k_i$. Thus every $k$ boils down to an ODE in $\tau$ where $k$ is a parameter, and each ODE is well-defined on its own. I have a custom type that wraps multiple ODE problems/solutions to take care of the interpolation in $k$ etc.

There are a lot of assumptions that have to be true for this to work out. The problem isn't that your case has a trivial mapping to many ODEs, it's that many other cases do not. You could have a term x(t, k+1), Differential(k)(x(t,k)), Integral(k)(x(t,k)), a non-local function of x(t) (such as average over k of x(t,k)), etc. You generally need some PDE discretization to turn the PDE into a solvable set of ODEs.

For this one, not only do you have that the equations are independent across k, but you also have a continuous and differentiable dependence in k. Thus while in general x(t,k) and x(t,k+1) can just be completely different, you also have the property that x(t,k+epsilon) is close to x(t,k). You're then using that property to say you don't need to solve at every k (which is computationally impossible, there are infinitely many), but you can instead just discretize some dk and solve x(t,k+dk*i) for i ODEs. But see that's a smoothness/regularity assumption in k which you'd need to prove: you cannot do this for general equations.

But that doesn't mean MTK doesn't cover this. You'd just phrase it as a PDESystem and then we'd need to add a discretize dispatch that is able to create the EnsembleProblem of ODEs to solve. This would be a reasonable discretizer to add as a default and is a feature request that should get an issue.

If you do it like that, then you'd also get the PDE plotting, so you could plot slices of x(t,:) and all sorts of other niceties built around the PDE form.

[hersle]

Thanks, this is very insightful. I see what you mean that MTK needs to cater for the general case and not the special case.

If you do it like that, then you'd also get the PDE plotting, so you could plot slices of x(t,:) and all sorts of other niceties built around the PDE form.

This is some functionality I am already making for my custom solution type, so tapping into this down the line would be very nice.

This would be a reasonable discretizer to add as a default and is a feature request that should get an issue.

I will close this PR for now, and we can resume this elsewhere. I do appreciate the discussion and help, as always. Thank you.