# How to perform integration processes NDSOlve and show list of random variables used in this process?

I have an ODE system which solves 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 r[n_]). In addition, the system has a other random variable (\[Xi][j_, k_] different of r[n_]). The code of the system is shown below:

nmax = 5;(*max number of equations*)
T = 70;(*Max Time*)
\[Xi][j_, k_] := RandomReal[];(*random variable*)
r[n_] := RandomInteger[{1,
n}];(*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] - (Sum[\[Xi][j, k] Subscript[x, k][t] Boole[k != j], {k,
n}])), {j, n}];
ics = {Subscript[x, 1][0] == 0.7};(*initial ICs*)
(*main loop*)
Do[(*solve for n species*)
sol[n] = NDSolve[{eqns, ics}, vars, {t, 0, T}][[1]];
(*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, r[n]][T] /. sol[n]]}];, {n, nmax}]


I try to do two things:

i) Do 1000 different process as described above;

I thought about adding a Do as:

Do[Do[(*solve for n species*)
sol[n] = NDSolve[{eqns, ics}, vars, {t, 0, T}][[1]];
(*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, r[n]][T] /. sol[n]]}];, {n, nmax}],{pp,1000}]


But I was not successful.

ii) Obtain a list with all the values of the variable \[Xi][j, k] drawn in the 1000 operations.

Can someone help me?

• "I was not successful" -- what went wrong? Mar 8, 2017 at 20:37
• @ChrisK The errors: NDSolve::ndode: Input is not an ordinary differential equation. ReplaceAll::reps: {(Subscript[x, 1]^[Prime])[t]==(1. -Subscript[x, 1][t]) Subscript[x, 1][t]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. ReplaceAll::reps: {Subscript[x, 1][0]==0.223725,Subscript[x, 2][0]==1.39449*10^-8,Subscript[x, 3][0]==0.000173822,Subscript[x, 4][0]==0.715508,Subscript[x, 5][0]==0.810497,Subscript[x, 6][0]==0.00223725} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.
– SAC
Mar 8, 2017 at 22:19

1. Your attempt in (i) fails because you forgot to initialize ics for every round of repeation.
2. Generating the list of ξ[j, k] before the summation will make the collection of it easier.
3. It's better to define vars and eqns as functions of n.

The following is my solution:

nmax = 5; T = 70; repeat = 1000;
r[n_] := RandomInteger[{1, n}];
vars[n_] := Table[x[j], {j, n}];
eqns[n_, ξ_] :=
Table[x[j]'[t] ==
x[j][t] (1 - x[j][t] - Sum[ξ[[j, k]] x[k][t] Boole[k != j], {k, n}]), {j, n}];

ans = Table[Module[{ics = {x[1][0] == 0.7}},
Table[
Module[{xi = RandomReal[1, {n, n}], sol, icvalue},
sol = NDSolveValue[{eqns[n, xi], ics}, vars[n], {t, 0, T}];
icvalue = Through@sol[T];
ics = {Thread[Table[x[j][0], {j, n}] == icvalue],
x[n + 1][0] == 0.01 icvalue[[r@n]]};
{sol, xi}], {n, nmax}]],
{repeat}];

• Thank you very much for your suggestion @xzczd
– SAC
Mar 10, 2017 at 16:00