# Random number in system of equations, solving using NDSolve & WhenEvent

I need to analyse the effect of random forcing on the system of coupled equations. For example, I have shown the equations below.

$$\begin{pmatrix} x1''\\x2''\\x3''\\x4''\\x5'' \end{pmatrix} + M_{(5 \times 5)} \begin{pmatrix} x1\\x2\\x3\\x4\\x5 \end{pmatrix} + V2_{(5 \times 1)} xb = \begin{Bmatrix} f1\\ f2\\ f3\\ f4\\ f5 \end{Bmatrix}, \\ \\ \text{and } \;\;\; xb=f(x5,t)$$

Here, F={f1,f2,f3,f4,f5} is the random force vector. I have used WhenEvent which induces random number at ceratain dt interval for random force. For two simplified coupled equations, I have mentioned the code below according to this answer.

tf=100; dt=0.001;
sol2=NDSolve[{x''[t]+x'[t]+2*x[t]-y[t]==r1[t],y''[t]+y'[t]+2*y[t]-x[t]==r2[t],WhenEvent[Mod[t,dt]==0,
{r1[t],r2[t]}->RandomReal[{-0.01,0.01},2]],{x[0]==y[0]==x'[0]==y'[0]==r1[0]==r2[0]==0}},{x,y},{t,0,tf},
DiscreteVariables->{r1,r2},Method->{"StiffnessSwitching",Method->{"ExplicitRungeKutta",Automatic}},
MaxSteps->Infinity,AccuracyGoal->7,PrecisionGoal->7][[1]];


Above code is working, but, I need to carryout the analysis for a large system of equations, say more than 20. So, for that I have written down the following code.

tf = 10; dt = 0.01;  N1 = 10;  \Omega=2;
Y1[t_] := Table[Subscript[x, i][t], {i, 2, N1}]; F[t_] := Table[Subscript[f, i][t], {i, 2, N1}];
a = ConstantArray[0, N1 - 1]; a[[N1-1]]=1;
xb[t] = Subscript[x, N1][t] \Cos[\Omega t];
eqns = Thread[D[Y1[t], t, t] + D[Y1[t], t] + Y1[t] ==a*xb[t]+ F[t]];
X0 = ConstantArray[0, N1 - 1]; Xd0 = ConstantArray[1, N1 - 1];
sol = NDSolve[{eqns, WhenEvent[Mod[t, dt] == 0, F[t] -> RandomReal[{-0.01, 0.01}, 2]],
Y1[0] == F[0] == X0, Y1'[0] == Xd0}, {Y1[t], F[t]}, {t, 0, tf}, DiscreteVariables -> F,
Method -> {"StiffnessSwitching",  Method -> {"ExplicitRungeKutta", Automatic}},
MaxSteps -> Infinity, AccuracyGoal -> 7, PrecisionGoal -> 7][[1]]

Plot[Evaluate[{ Subscript[x, N1][t]} /. First@sol], {t, 0, tf},
PlotRange -> Full]


It is working well in the absence of random number generator and WhenEvent approach. While consideration of whenevent for random number generation, it gives the following error.

NDSolve::underdet: There are more dependent variables, {Subscript[f, 2][t],Subscript[f, 3][t],Subscript[x, 2][t],Subscript[x, 3][t],{Subscript[f, 2][t],Subscript[f, 3][t]}}, than equations, so the system is underdetermined.

Can someone please help me, to solve the issue?

• The random force is the white noise, which I have managed with RandomReal function, as I believe that adapting standard stochastic differential equation process with the nonlinear system of equations will be very complex.
– Hari
Aug 9, 2020 at 14:24
• Can you include the error you mention? It seems this is missing from your question. Aug 10, 2020 at 9:08
• It is added. Thanks
– Hari
Aug 10, 2020 at 11:18
• Do your equations interact with each other or are they independent? Aug 12, 2020 at 17:11
• Here, I have shown the simplified form of system of equations. Actually the coefficient matrices are not diagonal. So, it is a coupled system of equation. a is also not a zero vector.
– Hari
Aug 13, 2020 at 6:01

# Updated to Accommodate Change in $$a$$ Definition

The OP changed the definition of the symbol, $$a$$ from my original answer. The code is essentially the same as before:

tf = 10; dt = 0.1; N1 = 10; Ω = 2;
a = ConstantArray[0, N1 - 1]; a[[N1 - 1]] = 1;
zeros = ConstantArray[0, N1 - 1];
Y1[t_] := Table[Subscript[x, i][t], {i, 2, N1}];
F[t_] := Table[Subscript[f, i][t], {i, 2, N1}];
dvbls = Head[#] & /@ F[t];
xb[t] = Subscript[x, N1][t] Cos[Ω t];
eqns = Thread[D[Y1[t], t, t] + D[Y1[t], t] + Y1[t] == a*xb[t] + F[t]];
ics = Thread[Y1[0] == zeros]~Join~Thread[Y1'[0] == zeros + 1];
istates = Thread[F[0] == zeros];
events = WhenEvent[Mod[t, dt] == 0,
Subscript[f, #][t] -> RandomReal[{-0.2, 0.2}]] & /@ Range[2, N1];
(* Join System of Equations to be Solved *)
system = Join[eqns, ics, istates, events];
sol = NDSolve[system, Join[Y1[t], F[t]], {t, 0, tf},
DiscreteVariables -> dvbls,
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}},
MaxSteps -> Infinity, AccuracyGoal -> 7, PrecisionGoal -> 7];
Plot[Evaluate[Table[Subscript[x, i][t], {i, 2, N1}] /. sol], {t, 0,
tf}, PlotRange -> Full]
Plot[Evaluate[Table[Subscript[f, i][t], {i, 2, N1}] /. sol], {t, 0,
tf}, PlotRange -> Full]


# Original Answer

I am not sure if your equations sense (e.g., in your case, the Cos term cancels because $$a$$ is a zero array) or just a simplified version, but you will need to Thread your initial conditions and states like you did for eqns. Also, you will need to create a list of events. Below is an example (I adjusted some parameters for visualization):

tf = 10; dt = 0.1; N1 = 10; Ω = 2;
a = ConstantArray[0, N1 - 1];
zeros = ConstantArray[0, N1 - 1];
Y1[t_] := Table[Subscript[x, i][t], {i, 2, N1}];
F[t_] := Table[Subscript[f, i][t], {i, 2, N1}];
dvbls = Head[#] & /@ F[t];
xb[t] = Subscript[x, N1][t] Cos[Ω t];
eqns = Thread[D[Y1[t], t, t] + D[Y1[t], t] + Y1[t] == a*xb[t] + F[t]];
ics = Thread[Y1[0] == zeros]~Join~Thread[Y1'[0] == zeros + 1];
istates = Thread[F[0] == zeros];
events = WhenEvent[Mod[t, dt] == 0,
Subscript[f, #][t] -> RandomReal[{-0.2, 0.2}]] & /@ Range[2, N1];
(* Join System of Equations to be Solved *)
system = Join[eqns, ics, istates, events];
sol = NDSolve[system, Join[Y1[t], F[t]], {t, 0, tf},
DiscreteVariables -> dvbls,
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}},
MaxSteps -> Infinity, AccuracyGoal -> 7, PrecisionGoal -> 7];
Plot[Evaluate[Table[Subscript[x, i][t], {i, 2, N1}] /. sol], {t, 0,
tf}, PlotRange -> Full]
Plot[Evaluate[Table[Subscript[f, i][t], {i, 2, N1}] /. sol], {t, 0,
tf}, PlotRange -> Full]


• Yes it's working. Very useful concept of joining the variables. Before updating also the code was perfect. Thank you very much
– Hari
Aug 15, 2020 at 8:23
• @Hari You are welcome. I'm glad that I could help. Aug 15, 2020 at 13:30