# How to solve an ODE system that periodically increases in size

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?

• I'm sorry but your question is unclear. What does it mean "ODE system that increases in size according to the rules"? and "This code solve the system of according with interpolation structure {{x1},{x1,x2},{x1,x2,x3},{x1,x2,x3,x4,},{x1,x2,x3,x4,x5}}"? Could you explain what problem you're trying to solve? Feb 23, 2017 at 20:34
• To solve the above system 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}] sol0 = Table[s = NDSolve[eqns[[l]], vars[[l]], {t, 0, T}], {l, n}] I forgot to add the NDSOlve part
– SAC
Feb 23, 2017 at 21:31
• @ChrisK The result of interpolation NDSolve provides the results of the form {{x1},{x1,x2},{x1,x2,x3},{x1,x2,x3,x4,},{x1,x2,x3,x4,x5}}. The initial conditions follows the same structure. I need to give the initial condition for each j by taking into account the last value given to j previous. For example, for j=1 Subscript [x, 1] [0] = 0.7. For j = 2, the initial conditions Subscript [x, 1] [0] = Subscript [x, 1] [T], Subscript [x, 2] [0]= Subscript [x, 1] [T] , Where Subscript [x, 1] [T] is the last value of the interpolation obtained at j = 1.
– SAC
Feb 23, 2017 at 21:45
• I think verbally describing the problem you want to solve will be clearer than using Mathematica code. Feb 23, 2017 at 21:57
• If you share the mathematical form of the problem then it might help?
– zhk
Feb 24, 2017 at 13:15

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


• Thanks for the answer. If possible, please include a phase portrait?
– zhk
Feb 25, 2017 at 18:54
• @MMM The dimension of the system keeps increasing, so that'd be tough without ParametricPlot5D! :) :) Feb 25, 2017 at 18:56
• As a test case, restrict it to 3 species or even 2.
– zhk
Feb 25, 2017 at 18:59
• @MMM OK see edit above Feb 25, 2017 at 19:07
• I did. Thankyou
– zhk
Feb 25, 2017 at 19:08

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}];

• Hi @Anjan this idea is exactly what I'm trying to solve. I tried to use your code suggestion, however the following errors appeared and I could not fix them. Do you know what is happening?
– SAC
Feb 24, 2017 at 12:16
• FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option. >> NDSolve::berr: There are significant errors {-1.,0.} in the boundary value residuals. Returning the best solution found. >>
– SAC
Feb 24, 2017 at 12:25
• FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option. >> NDSolve::berr: There are significant errors {-1.,0.} in the boundary value residuals. Returning the best solution found. >> FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option. >>
– SAC
Feb 24, 2017 at 12:26
• General::stop: Further output of FindRoot::sszero will be suppressed during this calculation. >> NDSolve::berr: There are significant errors {-1.,0.} in the boundary value residuals. Returning the best solution found. >> General::stop: Further output of NDSolve::berr will be suppressed during this calculation. >>
– SAC
Feb 24, 2017 at 12:27
• At what value of i is it not working? Feb 24, 2017 at 12:53