如何在Mathematica中模拟多个点电荷（球轴承）之间的斥力？

Question

我正在尝试在Mathematica中编写一个程序，模拟带有同样电荷（它们相互排斥）的带电滚珠扩散的方式。到目前为止，我的程序可以防止滚珠从屏幕上离开，并计算它们击中盒子边缘的次数。目前我让滚珠随机移动在盒子里，但我需要知道如何使它们相互排斥。以下是我的代码：

Manipulate[
 (*If the number of points has been reduced, discard  points*)
 If[ballcount < Length[contents], 
   contents = Take[contents, ballcount]];

 (*If the number of points has been increased, generate some random points*)
 If[ballcount > Length[contents], 
  contents = 
   Join[contents, 
    Table[{RandomReal[{-size, size}, {2}], {Cos[#], Sin[#]} &[
       RandomReal[{0, 2 \[Pi]}]]}, {ballcount - Length[contents]}]]];

 Grid[{{Graphics[{PointSize[0.02],

  (*Draw the container*)
  Line[size {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}], 
  Blend[{Blue, Red}, charge/0.3],
  Point[

   (*Start the main dynamic actions*)
   Dynamic[

    (*Reset the collision counter*)
    collision = 0;

    (*Check for mouse interaction and add points if there has been one*)
    Refresh[
     If[pt =!= lastpt, If[ballcount =!= 50, ballcount++]; 
      AppendTo[
       contents, {pt, {Cos[#], Sin[#]} &[
         RandomReal[{0, 2 \[Pi]}]]}]; lastpt = pt], 
     TrackedSymbols -> {pt}];

    (*Update the position of the points using their velocity values*)
    contents = Map[{#[[1]] + #[[2]] charge, #[[2]]} &, contents];

    (*Check for and fix points that have exceeded the box in Y
      direction, incrementing the collision counter for each one*)
    contents = Map[
      If[Abs[#[[1, 2]]] > size, 
        collision++; {{#[[1, 1]], 
          2 size Sign[#[[1, 2]]] - #[[1, 2]]}, {1, -1} #[[
           2]]}, #] &,
      contents];


    (*Check for and fix points that have exceeded the box in X 
      direction, incrementing the collision counter for each one*)
    contents = Map[
      If[Abs[#[[1, 1]]] > size, 
        collision++; {{2 size Sign[#[[1, 1]]] - #[[1, 1]], #[[1, 
           2]]}, {-1, 1} #[[2]]}, #] &,
      contents];

    hits = Take[PadLeft[Append[hits, collision/size], 200], 200];
    Map[First, contents]]]},
 PlotRange -> {{-1.01, 1.01}, {-1.01, 1.01}}, 
 ImageSize -> {250, 250}],

(*Show the hits*)
Dynamic@Show
  [
   ListPlot
   [
   Take[MovingAverage[hits, smooth], -100
    ]
   ,
   Joined -> True, ImageSize -> {250, 250}, AspectRatio -> 1, 
   PlotLabel -> "number of hits", AxesLabel -> {"time", "hits"}, 
   PlotRange -> {0, Max[Max[hits], 1]}], Graphics[]
  ]
}}
  ]
 ,
 {{pt, {0, 1}}, {-1, -1}, {1, 1}, Locator, Appearance -> None},
 {{ballcount, 5, "number of ball bearings"}, 1, 50, 1},
 {{charge, 0.05, "charge"}, 0.002, 0.3},
 {smooth, 1, ControlType -> None, Appearance -> None},
 {size, 1, ControlType -> None, Appearance -> None},
 {hits, {{}}, ControlType -> None},
 {contents, {{}}, ControlType -> None},
 {lastpt, {{0, 0}}, ControlType -> None}
 ]

- Fermat11235

---

你可能需要使用 NDSolve 处理运动方程。通过在每一步中使用Nearest 来近似处理一个近场，并仅使用距离给定粒子一定距离内的粒子来为其力场做出贡献。如果你只有几十个粒子，这可能并不重要。 - Daniel Lichtblau

Answer 1

你的模拟需要一个“碰撞检测算法”。这些算法的领域非常广泛，因为它与电脑游戏一样古老（如乒乓）且不可能在此完全回答。

由于你将充电球每个时间步骤都推进，所以你现在的模拟非常基础，并使它们从一个位置“跳”到另一个位置。如果运动像恒定速度和零加速度一样简单，则您可以通过简单地将时间放入公式中计算出所有位置，知道 移动的确切方程式。当球从墙上弹回时，它会得到新的方程。

有了这个，您可以预测两个球何时碰撞。您只需解决两个球是否同时具有相同位置即可。这称为先验检测。当您采用现在的模拟时，您必须在每个时间步检查是否有两个球如此接近，以至于它们可能相撞。

问题是，您的模拟速度并不是无限高的，而且球越快，模拟中的跳跃就越大。然后，两个球不太可能超过彼此，您错过了碰撞。

考虑到这一点，您可以先阅读维基百科文章了解该主题的概述。接下来，您可以阅读一些科学文章或查看裂缝是如何工作的。例如，花栗鼠物理引擎是一个令人惊叹的2D物理引擎。为确保此类内容正常工作，我相信他们必须对其碰撞检测进行了很多思考。

- halirutan

Answer 2

我无法完成检测部分，因为那是一个独立的项目。但现在交互速度更快了，你有一些不必要的动态效果（如果你看到单元格一直很忙碌，这意味着某个动态效果正在错误地运行）。

我还添加了一个停止/开始按钮。我可能没有完全理解你的所有操作，但足够让我进行所做的更改。你也在使用AppendTo。你应该尝试预先分配内容，并使用Part[]来访问它，这样会更快，因为你似乎知道允许的最大点数？

我喜欢更加分散代码，这有助于我更好地看到逻辑。以下是屏幕截图和更新版本代码。希望你发现它更快了。

请参见下面的代码，在更新（1）中。

(*updated version 12/30/11 9:40 AM*)
Manipulate[(*If the number of points has been reduced,discard points*)
\


 Module[{tbl, rand, npt, ballsToAdd},

  If[running,
   (
    tick += $MachineEpsilon;
    If[ballcount < Length[contents], 
     contents = Take[contents, ballcount]];

    (*If the number of points has been increased,
    generate some random points*)

    If[ballcount > Length[contents],
     (
      ballsToAdd = ballcount - Length[contents];
      tbl = 
       Table[{RandomReal[{-size, size}, {2}], {Cos[#], Sin[#]} &[
          RandomReal[{0, 2 \[Pi]}]]}, {ballsToAdd}];
      contents = Join[contents, tbl]
      )
     ];

    image = Grid[{
       {LocatorPane[Dynamic[pt], Graphics[{

           PointSize[0.02],(*Draw the container*)
           Line[size {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}],
           Blend[{Blue, Red}, charge/0.3],

           Point[(*Start the main dynamic actions*)

            (*Reset the collision counter*)
            collision = 0;

            (*Check for mouse interaction and add points if there has \
been one*)
            If[EuclideanDistance[pt, lastpt] > 0.001, (*adjust*)
             (
              If[ballcount < MAXPOINTS,
               ballcount++
               ];

              rand = RandomReal[{0, 2 \[Pi]}];
              npt = {Cos[rand], Sin[rand]};
              AppendTo[contents, {pt, npt}  ];
              lastpt = pt
              )
             ];

            (*Update the position of the points using their velocity \
values*)

            contents = 
             Map[{#[[1]] + #[[2]] charge, #[[2]]} &, contents];

            (*Check for and fix points that have exceeded the box in \
Y direction,incrementing the collision counter for each one*)

            contents = Map[
              If[Abs[#[[1, 2]]] > size,
                (
                 collision++;
                 {{#[[1, 1]], 
                   2 size Sign[#[[1, 2]]] - #[[1, 2]]}, {1, -1} #[[2]]}
                 ),
                (
                 #
                 )
                ] &, contents
              ];

            (*Check for and fix points that have exceeded the box in \
X direction,
            incrementing the collision counter for each one*)


            contents = 
             Map[If[Abs[#[[1, 1]]] > size, 
                collision++; {{2 size Sign[#[[1, 1]]] - #[[1, 1]], #[[
                   1, 2]]}, {-1, 1} #[[2]]}, #] &, contents
              ];


            hits = Take[PadLeft[Append[hits, collision/size], 200], 
              200];
            Map[First, contents]
            ]
           },
          PlotRange -> {{-1.01, 1.01}, {-1.01, 1.01}}, 
          ImageSize -> {250, 250}
          ], Appearance -> None
         ],(*Show the hits*)

        Show[ListPlot[Take[MovingAverage[hits, smooth], -100],
          Joined -> True, ImageSize -> {250, 250}, AspectRatio -> 1,
          PlotLabel -> "number of hits", AxesLabel -> {"time", "hits"},
          PlotRange -> {0, Max[Max[hits], 1]}
          ]
         ]
        }
       }
      ]
    )
   ];

  image
  ],

 {{MAXPOINTS, 50}, None},
 {pt, {{0, 1}}, None},
 {{ballcount, 5, "number of ball bearings"}, 1, MAXPOINTS, 1, 
  Appearance -> "Labeled", ImageSize -> Small},
 {{charge, 0.05, "charge"}, 0.002, 0.3, Appearance -> "Labeled", 
  ImageSize -> Small},
 Row[{Button["START", {running = True; tick += $MachineEpsilon}], 
   Button["STOP", running = False]}],
 {{tick, 0}, None},
 {smooth, 1, None},
 {size, 1, None},
 {hits, {{}}, None},
 {{contents, {}}, None},
 {lastpt, {{0, 0}}, None},
 {{collision, 0}, None},
 {image, None},
 {{running, True}, None},
 TrackedSymbols ->     { tick},
 ContinuousAction -> False,
 SynchronousUpdating -> True

 ]

- Nasser

---

你的代码中的Point构造在我的Mathematica版本上会产生错误。我不知道是什么原因导致了这个错误，但是将主要动态代码放在Point语句之外（例如执行类似于<code for updating points>;Point[Map[First, contents]]而不是Point[<code for updating points>;Map[First, contents]]）似乎可以解决它。 - Heike

@Heike，谢谢你告诉我这个问题。但这很奇怪，因为它在我的V8.04上运行正常，可能是当我将代码复制并粘贴到SO时，我破坏了一些东西。我记得在复制后微调了一下，使它看起来更好。也许我弄坏了它。现在我要做的是粘贴我所拥有的内容，这是之前没有更改过的，并且这一次，在我粘贴后不要触摸任何东西，如果你可以尝试新的版本，如果对你有用，那么就是粘贴错误，然后我可以删除早期版本。谢谢。 - Nasser

顺便提一下，从笔记本复制 M 代码到 SO，然后再将其复制回 M 的问题在于它会被搞乱并以与笔记本中不同的方式重新格式化。我认为下次从/到 SO/笔记本进行复制时应该使用 CODE 单元格，因为笔记本中的 CODE 单元格不应该有这种重新格式化的问题。 - Nasser

是的，我刚刚验证了一下，是因为我在复制代码后在SO中编辑时犯了一个愚蠢的错误。我不小心多加了一个“，”。我真是太傻了。上面更新（1）中的代码是正确的。我现在会删除之前的副本。（下次我会学习不要在粘贴代码后触碰SO缓冲区中的代码）。谢谢。 - Nasser