28

Please help me to simulate the movement of a particle inside a region with elliptical walls such that particle is reflected from the walls and continues to move.

A friend was able write code to simulate a particle represented by a Disk bouncing around inside a square, but we can't do it for an ellipse.

x = 0.5;
y = 0.5;

vx = 1;
vy = Pi/2;

step = 0.01;
radius = 0.05;

Animate[
  x = x + vx*step;
  y = y + vy*step;
  If[Abs[x - 1] <= radius || Abs[x] <= radius , vx = -vx];
  If[Abs[y - 1] <= radius || Abs[y] <= radius, vy = -vy];
  Graphics[{
    Cyan, Rectangle[{0, 0}, {1, 1}],
    Gray, Disk[{x, y}, radius],
    Point[{0.0, 0.0}], Point[{1.0, 1.0}]
  }],
  {t, 0, Infinity}
]
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Spizhen
  • 283
  • 3
  • 5

3 Answers3

47

Edit V10!

This is simple example what we can now do in real time!

R = RegionUnion @@ Table[Disk[{Cos[i], Sin[i]}, .4], {i, 0, 2 Pi, Pi/6.}];
R2 = RegionBoundary@DiscretizeRegion@R;

go[] := (While[r > .105, x += v; r = RegionDistance[R2, x]; Pause[.01]]; bounce[];)


bounce[] := With[{normal = Normalize[x - RegionNearest[R2, x]]},
  If[break, Abort[]];
  v = .01 Normalize[v - 2 v.normal normal];
  x = x + v;
  r = RegionDistance[R2, x]; go[]
  ]


x = {1, 0.};
pos = {x};
break = False;
v = .01 Normalize@{2, 1.};
r = RegionDistance[R2, x];


RegionPlot[R2, Epilog -> Dynamic@Disk[x, .1], AspectRatio -> Automatic]
Button["break at edge", break = True;]
go[]

enter image description here

This is an example, not perfect but nice enough to start.

V9

Unfortunately I don't have time to explain now. But take a look at wikipedia ellips site, tangent line part especially.

DynamicModule[{u = 0, t0, imp, v1, x0 = {0, .49}, v0 = {.5, -1.0}, t, a = 1, b = .5, 
              c, f1, f2},
 DynamicWrapper[
  Graphics[{ Thick, Scale[Circle[], {a, b}], AbsolutePointSize@7, Dynamic@Point[x0],
             Dashed, Thin, Dynamic@Line[{{x0, imp}, {imp, imp + Normalize@v1}, 
                                         {imp - normal, imp + normal}}]
           }, PlotRange -> 1.1, ImageSize -> 500, Frame -> True],
  Refresh[
    If[(#/a)^2 + (#2/b)^2 & @@ x0 < 1,
       x0 += v0;,
       x0 = imp + v1; v0 = v1; rec]
    , TrackedSymbols :> {}, UpdateInterval -> .001]]
  ,
  Initialization :> (
    c = Sqrt[a^2 - b^2]; v0 = Normalize[v0]/100; f1 = {-c, 0}; f2 = {c, 0};

rec := ({t0, imp} = {t, x0 + t v0
           } /. Quiet@NSolve[(#/a)^2 + (#2/b)^2 &amp; @@ (x0 + t v0) == 1. &amp;&amp; 
                              t &gt; 0, t, Reals][[1]];
normal = Normalize[Normalize[imp - f1] + Normalize[imp - f2]];

v1 = Normalize[v0 - 2 normal (v0.normal)]/100;(*bounce*));

rec)]

enter image description here

Community
  • 1
Kuba
  • 136,707
  • 13
  • 279
  • 740
  • @Kuba, what is imp ? – Spizhen Dec 18 '13 at 21:31
  • @Spizhen imp = x0 + t v0 where t is from NSolve. So it is next impact position. – Kuba Dec 18 '13 at 21:33
  • @Kuba, and what happens if the a space between the variables simply placed? Here, for example imp = x0 + t v0 I do not know about this. Maybe it's a feature version 9? – Spizhen Dec 18 '13 at 21:44
  • @Kuba, space between t and v0. what is it? – Spizhen Dec 19 '13 at 06:22
  • @Spizhen t is the parameter, time if you want, x0 and v0 are vectors, of position and velocity. imp = x0 + t v0 means that starting from x0 with velocity v0 the ball will hit the border after t time. – Kuba Dec 19 '13 at 06:57
  • Nice answer and nice update. I hadn't seen this before. +1 – Mr.Wizard Jan 27 '15 at 22:18
  • @Mr.Wizard Thanks! This gives so much fun! – Kuba Jan 27 '15 at 22:21
  • Yes, this is fun. How did you choose the example shape? It exhibits interesting behavior. – Mr.Wizard Jan 27 '15 at 22:27
  • @Mr.Wizard random positions were not interesting for me and this setup was the first I've thought about. I was not thinking about it much. – Kuba Jan 27 '15 at 22:31
  • @Kuba, sweet V10 code, I've just been generalising it to any number of particles, for which it continues to work great. However, after 510 bounces MMA hits a 1024 recursion limit and stops; I guess no-one watched the single ball above bounce 500 times! I guess this is because go calls bounce which in turns calls go again, etc. I've modified my n-ball code to remove the go function and simply have a single While loop encompassing the go code. I wondered if there are any reasons I'm missing which make your approach more favourable? – Quantum_Oli May 12 '16 at 17:03
  • @Quantum_Oli Thanks for kind words :) If I understood correctly, you moved go[] definition explicitly inside bounce? That's perfectly fine. I just like to keep code modularized so if go[] happens to be longer it will be more readable to keep it outside. You've motivated me to update this question, I have an idea how to slightly improve performance. – Kuba May 12 '16 at 17:54
  • @Kuba Sort of... There's still a separation between the go[] code and the bounce[] code. I've uploaded a screen cap if you'd like to see the modification for n balls. I've also used a RegionDistance function (df), I presume that's likely to factionally improve performance? Screen cap: http://imgur.com/Rjil0gv – Quantum_Oli May 12 '16 at 18:02
  • @Quantum_Oli creating RegionDistance function is great idea. I was going to change something more simple but I was distracted by life :P. So I'm checking distance to region at each step while it is safe to make Floor[currentDistance/dx] steps without that. Even though that closes distance is probably in different direction than v, it tells us that there is nothing close anyway. – Kuba May 12 '16 at 20:45
  • something broke in this code as of MMA v10.4.1. The RegionPlot returns only the empty rectangular frame, with the ball bouncing off invisible walls. Replacing that line with Show[R2, Graphics[Dynamic@{Disk[x, .1]}]] makes it work again (apart for some styling differences) – glS Aug 27 '16 at 12:06
24

My goal was quite ambitious. I wanted to create a way to let any rigid body bounce elastically against any other surface in MMA V9. To do this I use "masks" for the object and the environment. These masks are black and white images. White indicates that this is where the object/surface is, black is empty space. I can calculate the overlap between the object and the surface using Mathematica's image processing functions. Using the overlap I can calculate the normal of the surface. After that it's simple physics to change the velocity of the object accordingly. The code looks like this:

obj[mask_] := Graphics[{
   White, mask
   },
  PlotRange -> {{0, 500}, {0, 500}},
  ImageSize -> {500, 500},
  Background -> Black
  ]

forceVector[obj_, env_, center_] := N@Normalize[Plus @@ (center - # & /@ PixelValuePositions[ImageMultiply[obj, env], 1])]

step[{pt_, v_}] := Module[{f, nv},
  f = forceVector[obj[Disk[pt, 20]], ColorNegate@obj[Disk[{250, 250}, {100, 200}]], pt] /. (0. -> {0, 0});
  nv = If[v.f < 0, v - 2 v.f f, v];
  {pt + nv, nv}
  ]

pts = NestList[step, {{250, 250}, {1, 2}}, 1000];

frames = Graphics[{
     Black, Rectangle[{0, 0}, {500, 500}],
     White, Disk[{250, 250}, {100, 200}],
     Orange, Disk[#, 20]
     },
    PlotRange -> {{0, 500}, {0, 500}},
    ImageSize -> {500, 500}
    ] & /@ pts[[All, 1]];

ListAnimate[frames]

Here's a gif with a reduced number of frames:

bouncing disk

One can play with the velocity of the disk as well as the number of frames to get a longer path without as many calculations. This method is not very fast.

If you don't have the time/computing power to pre-compute the position list, you can still view the simulation using the code below. It will probably be very slow on many computers though (which is why I chose to pre-compute positions):

DynamicModule[{pt = {250, 250}, v = {6, 2}, f},
 Dynamic[
  f = forceVector[obj[Disk[pt, 20]], 
     ColorNegate@obj[Disk[{250, 250}, {100, 200}]], 
     pt] /. (0. -> {0, 0});
  If[v.f < 0, v = v - 2 v.f f];
  pt = pt + v;
  Graphics[{
    Black, Rectangle[{0, 0}, {500, 500}],
    White, Disk[{250, 250}, {100, 200}],
    Orange, Disk[pt, 20]
    },
   PlotRange -> {{0, 500}, {0, 500}},
   ImageSize -> {500, 500}
   ]
  ]
 ]
C. E.
  • 70,533
  • 6
  • 140
  • 264
  • Your method is very long build. Even if the decrease in the value of two times. Is this normal? How long should I wait? 10 minutes is not ready. P.S. Hardware: MacBook Air 13" (mid 2013) – Spizhen Dec 17 '13 at 20:13
  • @Spizhen Yeah, it's kind of slow because the way it detects collision is very expensive. I use an iMac and it certainly didn't take me ten minutes to execute that code but it may very well be the case on Macbook Air. – C. E. Dec 17 '13 at 20:23
  • whether there is a primitive example of solving tasks? Maybe a circle instead of an ellipse. Tomorrow morning need to show my teacher =( – Spizhen Dec 17 '13 at 20:32
  • 2
    @Spizhen whoah, aren't you ashamed writing such comments? have you tried anything? C. E., interesting approach ;) +1. – Kuba Dec 17 '13 at 20:35
  • @Kuba, I apologize, but I do not know what to do with it. Sorry guys – Spizhen Dec 17 '13 at 20:42
  • @Spizhen we are not offended, I just think that, since it is a homework, you should try to do something by yourself too – Kuba Dec 17 '13 at 20:45
  • @Kuba Sure. In any case, thank you. at the moment the problem is that the code is not compiled for a long time. That's why there is little panic =) – Spizhen Dec 17 '13 at 20:51
  • 1
    @Spizhen I added a version that you can show without having to do anything in advance. You still need the function definition for forceVector and obj from earlier. – C. E. Dec 18 '13 at 00:37
  • there is a semicolon missing after frames = ... – grbl Dec 18 '13 at 09:25
  • @grbl Thanks, updated. – C. E. Dec 18 '13 at 10:00
  • @Anon, Thank you. even as it turned out, in version 7 lacks a PixelValuePositions – Spizhen Dec 18 '13 at 14:51
  • @Spizhen Ah, I didn't know you were on version 7. It's new in version 9 :) – C. E. Dec 18 '13 at 16:50
14

  • @Kuba has provided an excellent solution. Here we follow his idea and use another approach like WhenEvent to get the particle tracing.

  • r[t] is the trace of the point.

Needs["OpenCascadeLink`"];
SeedRandom[1];
reflect[vector_, 
   normal_] = -(vector - 2 (vector - Projection[vector, normal])) // 
   Simplify;
shape = OpenCascadeShape[
   RegionUnion @@ 
    Table[Ball[{Cos[i], Sin[i], 0}, .4], {i, 0, 2 Pi, Pi/6.}]];
bm = OpenCascadeShapeSurfaceMeshToBoundaryMesh[shape, 
   "ShapeSurfaceMeshOptions" -> {"AngularDeflection" -> .1}];
R = BoundaryMeshRegion[bm];
R2 = RegionBoundary[R];
dist = RegionDistance[R2];
proj = RegionNearest[R2];
pt0 = RandomPoint[R, 1][[1]];
v0 = {1., 1., 0.};
d0 = 0.01*Norm[v0];
sol = NDSolveValue[{r''[t] == {0, 0, 0}, r[0] == pt0, r'[0] == v0, 
    WhenEvent[dist@r[t] <= d0, 
     r'[t] -> reflect[r'[t], r[t] - proj@r[t]]]}, r, {t, 0, 100}, 
   MaxStepSize -> 0.01];
ani = Animate[
  Show[Graphics3D[{{FaceForm[Opacity[.2], Darker@Cyan], EdgeForm[], 
      R2}}, Boxed -> False], 
   ParametricPlot3D[sol@t, {t, 0, c}, Mesh -> {{c}}, 
     MeshStyle -> {AbsolutePointSize[8], Red}, 
     Method -> {"BoundaryOffset" -> False}, 
     PlotStyle -> {Opacity[.9], White}, PlotPoints -> 400, 
     PerformanceGoal -> "Quality", PlotRange -> All] /. Line -> Arrow,
    ViewPoint -> {1, 1, .8}, Background -> Gray], {c, $MachineEpsilon,
    100}, AnimationRate -> 1, ControlPlacement -> Bottom]

enter image description here

  • Test another type of bounce and another obstacle surfaces.
SeedRandom[1]; 
reflect[vector_, 
  normal_] = -(vector - 2 (vector - Projection[vector, normal])) // 
  Simplify;
surf = ExampleData[{"Geometry3D", "UtahTeapot"}, "Region"];
cube = TransformedRegion[Cuboid @@ Transpose@RegionBounds[surf], 
     ScalingTransform[{1.2, 1.2, 1.2}, RegionCentroid[surf]]] // 
    DiscretizeRegion // RegionBoundary;
reg = RegionUnion[surf, cube];
dist = RegionDistance[reg];
proj = RegionNearest[reg];
pt0 = RegionCentroid[reg];
v0 = {1., 1., 0.2};
d0 = 0.01*Norm[v0];
sol = NDSolveValue[{r''[t] == {0, 0, -9.8}, r[0] == pt0, r'[0] == v0, 
   WhenEvent[dist@r[t] <= d0, 
    r'[t] -> reflect[r'[t], r[t] - proj@r[t]]]}, r, {t, 0, 20}, 
  MaxStepSize -> 0.001]
ani = Animate[
  Show[Graphics3D[{{FaceForm[Opacity[.2], Cyan], EdgeForm[], 
      surf}, {FaceForm[Opacity[.3], LightGray], EdgeForm[], cube}}, 
    Boxed -> False], 
   ParametricPlot3D[sol@t, {t, 0, c}, Mesh -> {{c}}, 
     MeshStyle -> {AbsolutePointSize[8], Red}, 
     Method -> {"BoundaryOffset" -> False}, 
     PlotStyle -> {Opacity[.9], White}, PlotPoints -> 400, 
     PerformanceGoal -> "Quality", PlotRange -> All] /. Line -> Arrow,
    Background -> LightGray], {c, $MachineEpsilon, 20}, 
  AnimationRate -> 1, ControlPlacement -> Bottom]

enter image description here

cvgmt
  • 72,231
  • 4
  • 75
  • 133