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.


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.


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.


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

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

  • $\begingroup$ What is n in the Sum and in the WhenEvent? $\endgroup$
    – Ivan
    May 28, 2015 at 0:58

1 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]

Mathematica graphics

  • $\begingroup$ 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)? $\endgroup$
    – march
    Aug 13, 2015 at 18:21
  • $\begingroup$ @march Yup. It's much faster (it takes out something like 75%) $\endgroup$ Aug 13, 2015 at 18:26
  • $\begingroup$ That's excellent. I definitely had a feeling that the way I was going about things would significantly slow things down. $\endgroup$
    – march
    Aug 13, 2015 at 18:28
  • $\begingroup$ how can I use this for two variable functions in a system of coupled pde? $\endgroup$ Apr 1, 2016 at 11:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.