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

Question

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} x1x2x3x4x5 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?

Tim Laska · Accepted Answer

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]

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

One Answer

Updated to Accommodate Change in $a$ Definition

Original Answer

Add your own answers!

Ask a Question