Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
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. –  rasher Mar 6 at 11:50

1 Answer 1

up vote 1 down vote accepted

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

share|improve this answer
    
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. –  BurgerAutomata Mar 6 at 13:32

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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