0
$\begingroup$

I am modelling 3D motion of particles inside a force-field using second-order coupled PDEs, which I solve from time t=t0 to t=t1 with initial conditions x0,y0,z0. I have a total number of "particles" n with subscript i and for each of these particles I model j collisions.

It is relatively straightforward to obtain the position vector r[t1]=={x[t1],y[t1],z[t1]} and the velocity vector r'[t1] for a single particle p_i by solving the PDEs once, but what I want to do is keep track of these outputs and use them as inputs in the subsequent calculation. That is to say, I want to then solve my PDEs from t=t1 to t=t2 using initial conditions r[t1] and r'[t1], and so on until I'm done . My code so far is:

f[nn_] := 
Module[{n = nn}, 
  temp = Quiet@
    ParallelTable[(NDSolve[{x''[t] == -x'[t] - x[t], 
         y''[t] == -y'[t] - y[t], z''[t] == -z'[t] - z[t], x[0] == 1, 
         x'[0] == 1, y[0] == 1, y'[0] == 1, z[0] == 1, 
         z'[0] == 1}, {x, y, z}, {t, tjs[[i, j]], tjs[[i, j + 1]]}]) //
       First, {i, n}, {j, Length[tjs[[i]]] - 1}]]
nn = 20;
et = 2.1 10^7;
p = 1;
pdft[t_, p_] = et p 10^-6 E^(-(et p 10^-6) t);
trndm[p_] := 
 RandomVariate[ProbabilityDistribution[pdft[t, p], {t, 0, Infinity}]]
tj := Prepend[(sum = 0;
   Reap[While[sum < 0.05, sum = Sow[sum + trndm[p]]]][[2, 1]]), 0]
tjs = Table[tj, nn];
Table[(tjs[[i, Length[tjs[[i]]]]] = 0.05), {i, nn}];
f[nn]

So, each time I run the module f[nn] it properly solves j collisions for each particle i, but it does so using the same initial conditions every time.

I've tried to provide a MWE so if something is missing or needs clarification please let me know.

Improve this question
$\endgroup$
5
  • $\begingroup$ Your code does not work as posted (NDSolve complains that the system is underdetermined). 1) copy and past the code you posted, see if you can fix it so it runs; 2) for now, do away with the complication of parallel execution and concentrated on getting the code to run. $\endgroup$
    – MarcoB
    Mar 6, 2019 at 22:15
  • $\begingroup$ I intended to provide the algorithm without going into details, but I'll add a version of the code that actually executes now. $\endgroup$
    – Milad Dakka
    Mar 7, 2019 at 1:45
  • $\begingroup$ @MarcoB, it is done! $\endgroup$
    – Milad Dakka
    Mar 7, 2019 at 1:52
  • $\begingroup$ @MiladDakka What do you want to get in the result? $\endgroup$
    – Alex Trounev
    Mar 7, 2019 at 11:56
  • $\begingroup$ @AlexTrounev, I have answered my own question! What I personally want is to perform a statistical analysis of the result of millions of collisions for a particular research question I am investigating. For this purpose, I needed to "follow" the trajectory of atoms for a certain period of time, a trajectory that includes several "collisions", which I modelled as velocity changes. See the answer below for more info. $\endgroup$
    – Milad Dakka
    Mar 19, 2019 at 11:02

1 Answer 1

Reset to default
1
$\begingroup$

I have found an answer to my own question using advice from a Mathematica guru and I had to update this question as a result.

So, it turns out the WhenEvent function in Mathematica is the most wonderful answer to this question, as it allows one to specify an event (in my case a collision occurring at pre-determined times) and what one wants to happen as a result (in my case, the atom experiences a velocity change due to the collision) neatly within the NDSolve function.

So my example code, which includes "collisions" at times t=0.01, t=0.05 and t=0.08, now looks something like this:

NDSolveValue[
 {
  x''[t] == -x'[t] - x[t], 
  y''[t] == -y'[t] - y[t], 
  z''[t] == -z'[t] - z[t], 
  x[0] == 1, x'[0] == 1, 
  y[0] == 1, y'[0] == 1, 
  z[0] == 1, z'[0] == 1, 
  WhenEvent[
   {t==0.01,t==0.05,t==0.08},
   {x'[t]->RandomReal[{0,1}],
   y'[t]->RandomReal[{0,1}],
   z'[t]->RandomReal[{0,1}]}
  ]
 },
{x,y,z},
{t,0,0.1}
]

Here is a very helpful link to Events and Discontinuities in DEs, with examples for many different cases that could be of use to anyone who comes across this kind of issue.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Are you sure this code describes a collision of molecules? $\endgroup$
    – Alex Trounev
    Mar 19, 2019 at 11:36
  • $\begingroup$ Hi Alex, this code isn't meant to represent real collisions, only the structure of the code that could do so. The forces in this case are not very realistic and the collisions are assumed to result in random velocities. But yes, the essential idea is useful for describing collisions (or events). $\endgroup$
    – Milad Dakka
    Mar 24, 2019 at 10:06
  • 1
    $\begingroup$ There are several posts on this topic, for example mathematica.stackexchange.com/questions/38687/… $\endgroup$
    – Alex Trounev
    Mar 24, 2019 at 11:22

Your Answer

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

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