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

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".

• What is n in the Sum and in the WhenEvent? – Ivan May 28 '15 at 0:58

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]


• Thanks! I'll play around with this (and check myself), but did you by chance check to see if it's faster than the original code (or perhaps it will allow larger n's, even if it doesn't run faster)? – march Aug 13 '15 at 18:21
• @march Yup. It's much faster (it takes out something like 75%) – Dr. belisarius Aug 13 '15 at 18:26
• That's excellent. I definitely had a feeling that the way I was going about things would significantly slow things down. – march Aug 13 '15 at 18:28
• how can I use this for two variable functions in a system of coupled pde? – Ahmad Sheikhzada Apr 1 '16 at 11:23