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.