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.

This Code produces a box filled with balls at different positions with random charges on them and then calculates their motion according to Coulombs law. rendering the electric field as well as the potential at any point as a density plot along with it, it works just fine, but it doesnt seem to do the greatest of jobs at confining the balls to the box (it works about half of the time, i reccon if forces on a ball at the boundaries gets to big it accelerates past the "boundary trigger" if that makes any sense. I might be totally off too

as always, any help is appreciated!

enter image description here

this is how the rendered plot should look like for a single instance in time

g = 1;
Tmax = 30;
n = 10;
pos = RandomReal[{-10, 10}, {n, 2}]
charge = RandomReal[{-10, 10}, n]

Graphics[
 Table[
  Disk[pos[[i]], 0.5]
  , {i, 1, n}]
 , PlotRange -> 10.5]

sol = NDSolve[
  Flatten@Table[
    {
     With[{i = i}, {
       (y[i])''[t] == -g + VectorField[x[i][t], y[i][t], t][[2]],
       (y[i])'[0] == 0,
       y[i][0] == pos[[i, 2]],
       (x[i])''[t] == VectorField[x[i][t], y[i][t], t][[1]],
       (x[i])'[0] == 0,
       x[i][0] == pos[[i, 1]],
       WhenEvent[y[i][t] == -10, y[i]'[t] -> -0.95*y[i]'[t]],
       WhenEvent[y[i][t] == 10, y[i]'[t] -> -0.95*y[i]'[t]],
       WhenEvent[x[i][t] == -10, x[i]'[t] -> -0.95*x[i]'[t]],
       WhenEvent[x[i][t] == 10, x[i]'[t] -> -0.95*x[i]'[t]]
       }
      ]
     }
    , {i, 1, n}]
  ,
  Flatten@Table[
    {
     x[i], y[i]
     }
    , {i, 1, n}]
  , {t, 0, Tmax}, PrecisionGoal -> 4]

Manipulate[Graphics[
  Table[
   With[{p = {x[i][t], y[i][t]} /. First[sol]},
    {ColorData["TemperatureMap", (charge[[i]] + 10)/20], Disk[p, 0.5]}
    ]
   , {i, 1, n}]
  , PlotRange -> 10.5, Background -> Black]
 , {t, 0, Tmax}]

Fancy version:

PotentialAnimated[u_, v_, t_] := 
 Block[{x, y}, 
  Sum[charge[[i]]/
    EuclideanDistance[{u, v}, {x[i][t], y[i][t]} /. sol[[1]]], {i, 1, 
    n}]]

VectorField[u_, v_, t_] := 
 Block[{x, y}, 
  Sum[charge[[i]]*({u, v} - {x[i][t], y[i][t]} /. sol[[1]])/
     EuclideanDistance[{u, v}, {x[i][t], y[i][t]} /. sol[[1]]], {i, 1,
     n}]]

Manipulate[
 Show[{DensityPlot[
    PotentialAnimated[x, y, t], {x, -10, 10}, {y, -10, 10}, 
    ColorFunction -> "SunsetColors"], Graphics[
    Table[
     With[{p = {x[i][t], y[i][t]} /. First[sol]},
      {ColorData["TemperatureMap", (charge[[i]] + 10)/20], 
       Disk[p, 0.5]}
      ]
     , {i, 1, n}]
    , PlotRange -> 10.5, Background -> Black], 
   VectorPlot[{VectorField[x, y, t][[1]], 
     VectorField[x, y, t][[2]]}, {x, -10, 10}, {y, -10, 10}]}], {t, 0,
   Tmax}]
share|improve this question
1  
1  
What I would do is integrating for a small time step, then do the detection of collision outside NDsolve not inside using events, then call NDSolve again for a small time step, and so on. see this demo demonstrations.wolfram.com/… –  Nasser Jul 7 at 17:28
add comment

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.