How to solve an ODE system that periodically increases in size

Question

I have an ODE system that increases in size according to the rules

n = 5;
T = 50;
nu = 0.05;
vars = Table[Subscript[x, j][t], {i, n}, {j, i}];
eqns = Table[{Subscript[x, j]'[t] == 
    Subscript[x, j][
      t] (1 - Subscript[x, j][t] - 
       nu (Sum[ Subscript[x, k][t] Boole[k != j], {k, i}]) ), 
   Subscript[x, j][0] == 
    If[j == 1 && i == 1, 0.7, 
     If[j == i, 
      0.01 Subscript[x, RandomInteger[{1, j - 1}]][t] /. t -> T, 
      Subscript[x, j][t] /. t -> T]]}, {i, n}, {j, i}]

I want solve the differential equation of n variables, with initial conditions defined using the previous differential equation solution of n-1 variables and with an initial condition for the last variable (which randomly depends on one of the previous variables).

I'm having trouble building the initial conditions, present in the code above

 Subscript[x, j][0] == 
        If[j == 1 && i == 1, 0.7, 
         If[j == i, 
          0.01 Subscript[x, RandomInteger[{1, j - 1}]][t] /. t -> T, 
          Subscript[x, j][t] /. t -> T]]

Where the first initial condition {x1[0]=0.7}

The next {x1[0]=x1[Last time in j past],x2[0]=[Last time in j past]}

. . .

{{x1[0]=[Last time in j past],x2[0]=x2[Last time in j past],x3[0]==x3[Last time in j past],x4[0]=x4[Last time in j past],x5=0.01xSubscript[RandomInteger[1,j-1]][Last time in j past]}

can anybody help me?

Answer 1

This looks like a Lotka-Volterra competition model with weak, symmetric interspecific competition. If I understand correctly, you want to build up a community by introducing new species every T=50 time steps, where the new species has initial density of 0.01 times one of the existing species.

Call me old-fashioned, but I think this is easiest to understand when handled in an iterative Do loop:

nmax = 5; (* max number of species *)
T = 50; (* period *)
nu = 0.05; (* interspecific competition coefficient *)

(* set up unknown vars and differential equations, for n species *)

vars := Table[Subscript[x, j], {j, n}];
eqns := Table[Subscript[x, j]'[t] == Subscript[x, j][t]
  (1 - Subscript[x, j][t]- nu (Sum[Subscript[x, k][t] Boole[k != j], {k, n}]))
, {j, n}];

(* initial ICs *)
ics = {Subscript[x, 1][0] == 0.7};

(* main loop *)
Do[
  (* solve for n species *)
  sol[n] = NDSolve[{eqns, ics}, vars, {t, 0, T}][[1]];
  (* plot dynamics *)
  Print[Plot[Evaluate[Table[Subscript[x, j][t], {j, n}] /. sol[n]], {t, 0, T}, PlotRange -> {0, All}]];

  (* set up ICs for n=n+1 species *)
  ics = Join[
    Table[Subscript[x, j][0] == Evaluate[Subscript[x, j][T] /. sol[n]], {j, n}],
    {Subscript[x, n + 1][0] == Evaluate[0.01 Subscript[x, RandomInteger[{1, n}]][T] /. sol[n]]}
   ];
, {n, nmax}]

Not sure how interesting this is, since I believe this symmetric LV system has a globally stable equilibrium at $x_i=1/(1-(n-1)nu)$ as long as $0<nu<1$, but you could use this as a basis for more interesting explorations of community assembly.

Edit

Here's a 3D phase portrait of the first three periods for @MMM:

Show[
  ParametricPlot3D[Evaluate[{Subscript[x, 1][t], 0, 0} /. sol[1]], {t, 0, T}],
  ParametricPlot3D[Evaluate[{Subscript[x, 1][t], Subscript[x, 2][t], 0} /. sol[2]], {t, 0, T}],
  ParametricPlot3D[Evaluate[{Subscript[x, 1][t], Subscript[x, 2][t], Subscript[x, 3][t]} /. sol[3]], {t, 0, T}],
  PlotRange -> {{0, 1}, {0, 1}, {0, 1}}
]

Answer 2

If I understand it correctly, I think the following is what you want:

Solve the differential equation of n variables, with initial conditions defined using the previous differential equation solution of n-1 variables and with an initial condition for the last variable (which randomly depends on one of the previous variables).

The following code uses iniCond which generates initial conditions for a given i, diffGen generates equations with the initial conditions, and finally solveDiffGen solves the differential equations.

T = 50;
nu = 0.05;
vars[i_] := Subscript[x, #][t] & /@ Range[i];
iniCond[1] = {Subscript[x, 1][0] == 0.7};
diffGen[i_] := 
  Join[Table[{Subscript[x, j]'[t] == 
       Subscript[x, j][
         t] (1 - Subscript[x, j][t] - 
          nu (Sum[Subscript[x, k][t] Boole[k != j], {k, i}]))}, {j, 
      i}], iniCond[i]] // Flatten;
iniCond[i_] := 
  Module[{prevSolConds, lastCond}, 
   prevSolConds = 
    NDSolve[diffGen[i - 1], 
       vars[i - 1], {t, 0, 
        T}] /. (x_ -> g_) :> {(x /. {t -> 0}) == (g /. {t -> T})} // 
     Flatten;
   lastCond = 
    Subscript[x, i][0] == 
      0.01 Subscript[x, RandomInteger[{1, i - 1}]][t] /. t -> T;
   {prevSolConds, lastCond} // Flatten];
solveDiffGen[i_] := NDSolve[diffGen[i], vars[i], {t, 0, T}];