2
$\begingroup$

I've been toying around with NDSolve for a while, and read through the website. By doing so I discovered that I could use it for vectors and arrays with much pleasure. So I wanted to write a simple simulator with multiple balls or particles interacting with Lennard Jones potentials. When I write this way I have no problem

partsim[p0_, v0_, mass_] := 
  Module[{rp},
  rp[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
  rp[t] /. NDSolve[
  Join[
   Table[
    {
     (*Initial COnditions*)
     rp[0][[i]] == p0[[i]],
     rp'[0][[i]]  == v0[[i]],

     (*Newton Equations*)  
     rp''[t][[i]] == (1/mass[[i]])*
       Sum[If[j != i, ljForce[i, j, rp[t]], 0], {j, 1, 
         nParticles}]
     (*Third newton law in the future*)

     }
    , {i, nParticles}](*End of Table*)  
   ](*End of Join*)
  , rp[t]
  , {t, 0, ftime}] // First
   ](*End of Function*);

The real problem arises when I'm trying to also add a WhenEvent for each of the particles present in the system. Since it's 3d I would like them to stay inside a bounding box and invert the sign of their velocity every time they touch the walls or at least come withing tolerance.

partsim[p0_, v0_, mass_] := 
  Module[{rp, vel},
  rp[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
  vel[t_] := Evaluate@Array[Unique[][t] &, {nParticles, 3}];
  rp[t] /. NDSolve[
  Join[
   Table[
    {
     (*Initial COnditions*)
     rp[0][[i]] == p0[[i]],
     rp'[0][[i]] == vel[0][[i]] == v0[[i]],
     (*Newton Equations*)
     rp'[t][[i]] == vel[t][[i]],        
     rp''[t][[i]] == (1/mass[[i]])*
       Sum[If[j != i, pairForce[i, j, rp[t]], 0], {j, 1, 
         nParticles}]


     }
    , {i, nParticles}](*End of Table*)


   (*Discrete Events*)
   (*Bouncing on walls*)

   , Table[{
    WhenEvent[
     {Norm[rp[t][[u]] - wallcoord] < .01  , 
      Norm[rp[t][[u]]] <  .01}
     , vel[t][[u]] -> -vel[t][[u]]
     ](*End of WhenEvent*)
    }, {u, nParticles}]

   ](*End of Join*)

  , rp[t]
  , {t, 0, ftime}, DiscreteVariables -> {vel}] // First
   ](*End of Function*);

The error that appears seems to be related to the fact that u remains unevaluated, as attempting to Table it outside of the program will reveal. I have added DiscreteVariables -> {vel} because I thought that since we are subjecting velocity to a sudden change it would be necessary. Is there a way of adding a list of working WhenEvent, one for each particle, to the system of equations? Many thanks

$\endgroup$
1
  • $\begingroup$ Please post a minimally working example, with everything required to use it - any data, parameters, arguments to use, etc.. As it stands, it is an exercise in mind-reading. $\endgroup$
    – ciao
    Mar 6, 2014 at 11:50

1 Answer 1

4
$\begingroup$

In lieu of a working example in the current OP to work with, yes, this can be done. An example of two bouncing balls with differing velocity retention on bounce:

eventtab = 
 Table[With[{n = n}, 
   WhenEvent[y[n][t] == 0, y[n]'[t] -> (-0.95 + n/100) y[n]'[t]]], {n,
    2}]

NDSolve[{y[1]''[t] == -9.81, y[1][0] == 5, y[1]'[0] == 0, 
  y[2]''[t] == -4, y[2][0] == 3, y[2]'[0] == 0, eventtab}, {y[1], 
  y[2]}, {t, 0, 10}]

Plot[Evaluate[{y[1][t], y[2][t]} /. %], {t, 0, 10}]

enter image description here

$\endgroup$
1
  • $\begingroup$ This was a very useful answer, I did not think about using With. I knew that WhenEvent had the Hold attribute, but I did not know I could bypass it. $\endgroup$ Mar 6, 2014 at 13:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.