Using a Mathematica index as a DiscreteVariable in NDSolve when solving a coupled set of ordinary differential equations

Question

Context

Since the explanation below of the problem to be solved is lengthy, let me preamble this by saying that I have code that works to solve the problem, but I don't know whether (1) it's optimized, (2) it uses best Mathematica practices. Specifically, (1) is my implementation of the sum using UnitStep functions (explained below) fast---i.e. is there a problem with how many evaluations are done when most of those evaluations end up not being used---and (2) is there a less kluge-y way to do this (not involving the sum over all the functions multiplied by UnitSteps that takes advantage of built-in, optimized Mathematica functions.

Question

I am solving a set of coupled ordinary differential equations in which the functions $f_i$ are indexed by a quantity $i$. There is a quantity $i^*$ that breaks the functions into two sets such that the $f_i$'s where $i\leq i^*$ satisfy one differential equation and the $f_i$'s where $i>i^*$ satisfy another. The right-hand sides of the equations for $i\leq i^*$ depend on all of the functions indexed by $i\leq i^*$. Finally, and most importantly, $i^*$ is itself a function of the $f_i$'s, so the differential equations are changing on the fly.

Example

As an example, let's consider the following prototype. For $i > i^*$, $f_i(t)$ exponentially decays. For $i \leq i^*$, $f_i(t)$ grows "exponentially" with a rate equal to $\sum_{j=1}^{i^*}f_j(t)$. Finally, $i^*(t) = \lfloor 1/f_1(t)\rfloor$, so that as $f_1(t)$ grows, $i^*$ decays. The differential equations are then given by \begin{align} \frac{df_i}{dt} &= \begin{cases} -f_i & i^* < i \leq N\\ f_i\sum_{j=1}^{i^*}f_j & i\leq i^* \end{cases}~, \\ i^* &= \left\lfloor \frac{1}{f_1}\right\rfloor~, \end{align} where $N$ is the number of functions to be solved for, and let's take as initial conditions \begin{align} f_{1\leq i \leq N}(0) &= \frac{1}{N}~,\\ i^*(0) &= N~. \end{align} Since $i^*(0) = N$, all of the functions will initially grow, and $i^*$ will decay until $1/f_1 = N-1$, at which point $f_N$ will exponentially decay from then on, etc.

The issue with implementing this set of equations with NDSolve is that the summation $\sum_{j=1}^{i^*}f_j$ cannot work. My initial thought was to use Sum with $i^*$ as the maximum index, and treat $i^*$ as a DiscreteVariable. This doesn't work. Following the working examples in this previous question along with a suggestion from a friend, I replaced the truncated sum with a sum over all the functions, with UnitStep functions that turn on and off the $f_i$ depending on how $i$ compares to the current value of $i^*$, and then treat $i^*$ as a DiscreteVariable. The problem is that if $i^*$ is relatively small, then that sum over all of the $f_i$ uses a lot of values that end up being multiplied by UnitStep functions that evaluates to zero, and hence it seems like there are a lot of unnecessary evaluations. Note that $N$ can be very large, on the order of at least hundreds.

Attempt

What follows is a simple example of working code that solves the prototype equation for arbitrary $N$.

numFunctions=3;
NDSolve[
 Join[
  Table[f[i]'[t] == Piecewise[{
    {-f[i][t], i > iS[t]}
    , {f[i][t] Sum[f[j][t] UnitStep[iS[t] - j], {j, numFunctions}] ,i <= iS[t]}
   }]
   , {i, numFunctions}]
  , Table[f[i][0] == 1/numFunctions, {i, numFunctions}]
  , {iS[0] == numFunctions}
  , {WhenEvent[Floor[1/f[1][t]] + 1 <= iS[t], iS[t] -> Floor[1/f[1][t]]]}
 ]
 , Join[Table[f[i], {i, numFunctions}], {iS}]
 , {t, 0, 10}
 , DiscreteVariables -> {iS}
]

Thanks in advance! Hopefully this is an "answerable question" rather than a "discussion-generating question".

Answer 1

Here you have a way without using neither UnitStep nor PieceWise that improves the performance by 75% wrt your code. It computes 200 functions in a very reasonable time for your toy example.

The main trick is to use a numeric (black box) function to be able to take Part[... ] inside it.

n = 200;
af = Array[f, n];
taf[t_] := Through[af[t]]
bb[i_?IntegerQ, lim_, f__] := f[[i]] If[i > lim, -1, Tr@f[[;; lim]]];

eqs = Join[Array[f[#]'[t] == bb[#, iS[t], taf@t] &, n],
           Thread[taf@0 == 1/n], {iS[0] == n},
           {WhenEvent[Floor[1/f[1][t]] + 1 <= iS[t], iS[t] -> Floor[1/f[1][t]]]}];
vars = Join[af, {iS}];

xsol = NDSolve[eqs, vars, {t, 0, 10}, DiscreteVariables -> {iS ∈ Integers}];

Plot[Array[f[#][t] &, n] /. xsol, {t, 0, 10}, PlotRange -> All, 
     AxesOrigin -> {0, 0}, Evaluated -> True]