Реализация базового метода Стёрмера-Верле / Хабр

Используем силу уравнений Ньютона и численных методов для моделирования динамики простых плоских мешей в реальном времени! В конце вы сможете моделировать падение ножниц ✂️ как на анимации

Зачем? Это наглядный, простой и интересный метод решения задачи многих тел с простым взаимодействием типа точка-палка-точка, что придется как раз для пятничного развлекательного программирования (recreational programming session!).

Среда разработки

Все примеры показаны в WLJS Notebook, так как автору лень изучать ipywidgets для создания простой интерактивщины, а Marimo слишком медленный.

Блокнот-версия доступна по ссылке

Набор точек

Создадим случайный набор точек:

pts0 = RandomReal[{-1,1}, {10, 2}];
pts1 = pts0; pts2 = pts0;

Graphics[{Rectangle[{-10,-10}, {10,10}], Red, Point[pts0 // Offload]}]

Чтобы понять, где будут эти точки в следующий момент времени, нам нужно решать как минимум, задачу кинематики, а так как в системе присутствуют внешние силы (гравитация), то и задачу динамики тоже. Метод Стёрмера–Верле — это лишь одна из численных схем для дискретизации уравнений движения Ньютона (II з.):

$x_{i+1} = 2x_i - x_{i-1} + a_i \delta t^2$

Чем она отличается от всех других? Для неё требуется хранить только текущую и две предыдущие позиции: $x_{i+1}$ , и $x_{i-1}$ . Как видно: скорость не нужна, что сыграет особую роль позже. Попробуйте выполнить следующий код:

With[{dt = 0.5},
  pts0  = 2 pts1 - pts2 + Table[{0,-1}, {Length[pts0]}] Power[dt,2];
  pts2 = pts1;
  pts1 = pts0;
]

Добавим случайные связи между точками (связи — это как твёрдые стержни):

bonds = RandomGraph[Length[pts0]{1, 2}, VertexCoordinates->pts0]

Метод Cтёрмера–Верле особенно удобен при наличии ограничений типа жёстких/мягких связей (как выше). Например, если хотим сохранить фиксированные расстояния между точками согласно графу для этого потребовалась бы нетривиальная работа со скоростью и координатой, но метод Верле оставляет нам только

Жесткость таких связей задаётся не через дифференциальные уравнения, а специальным алгоритмом, который аппроксимирует реальные потенциалы этих связей. На этот счёт есть целая наука, как не решать диффуры, а обойтись манипуляциями координат. Здесь же представлен наиболее наивный и простой вариант.

Пусть начальная длина "стержня" соединяющего две точки — $L = |\mathbf{d}|$ , где

$\mathbf{d} := \mathbf{x}_a - \mathbf{x}_b$

В какой-то момент при движении точек это расстояние будет отличным от того, что было в предыдущий момент времени. Тогда мы можем его нахально скорректировать

$\Delta \mathbf{d} := \mathbf{d} \cdot k \left( \frac{L}{|\mathbf{d}|} - 1 \right)$ $\begin{align*} \mathbf{x}_a &= \mathbf{x}_a + \frac{1}{2} \Delta \mathbf{d} \\ \mathbf{x}_b &= \mathbf{x}_b - \frac{1}{2} \Delta \mathbf{d} \end{align*}$

где k это параметр задающий жёсткость стержня (связи):

бесконечно жесткая палка
мягкая связь (как будто пружина aka закон Гука)

Реализация на чистых функциях:

getDir[bonds_, pts_] := Map[(pts[[#[[1]]]] - pts[[#[[2]]]]) &, bonds] 

applyBond[pairs_, k_:1.0, maxDelta_:0.1][p_] := Module[{pts = p},
  MapThread[Function[{edge, length}, With[{
    diff = pts[[edge[[1]]]] - pts[[edge[[2]]]]
  },
    With[{
      delta = Min[k ( ((length)/(Norm[diff]+0.001)) - 1), maxDelta] diff
    },
      pts[[edge[[1]]]] += delta/2.0;
      pts[[edge[[2]]]] -= delta/2.0;
    ]
  ] ], RandomSample[pairs] // Transpose];

  pts
]

Зададим функция для динамической визуализации наших точек и стержней

showGeometry[points_, bonds_, opts___] := Graphics[{
  Gray, Rectangle[5{-1,-1}, 5{1,1}], Blue,
  
  Table[
    With[{i = edge[[1]], j = edge[[2]]},
      Line[With[{p = points}, {p[[i]], p[[j]]}]] // Offload
    ]
  , {edge, bonds}],
  
  Red, Point[points // Offload]
}, opts]

SetAttributes[showGeometry, HoldFirst];

Нарисуем наши точки со стержнями

showGeometry[pts0, EdgeList @ bonds]

Теперь можно собрать основной цикл аним��ции:

pts0 = RandomReal[{-1,1}, {10, 2}];
pts1 = pts0; pts2 = pts0;

bonds = EdgeList[RandomGraph[Length[pts0]{1, 3}]];
lengths = Norm /@ getDir[bonds, pts0];
pairs = {bonds, lengths} // Transpose;

Do[With[{\[Delta]t = 0.05},
  pts0 = Clip[applyBond[pairs, 1.0][
    2 pts1 - pts2 + Table[{0,-1}, {i, Length[pts0]}] Power[\[Delta]t,2]
  ], {-5,5}];
  
  pts2 = pts1;
  pts1 = pts0;
  
  Pause[0.01];
], {i, 200}]

Процедурная генерация мешей

Даже без оптимизации (например, векторизации или OpenCLLink) можно попробовать более сложные конструкции и взаимодействие с мышью пользователя.

Обод

Создаём обод из точек:

rim = Join @@
  (Table[#{Sin[x], Cos[x]}, {x,-Pi, Pi, Pi/6}] &/@ {1, 0.9} // Transpose);
  
rbonds = Graph[
  Join[
    Table[i<->i+1, {i, 1, Length[rim]-1}],
    Table[i<->i+2, {i, 1, Length[rim]-2, 2 }],
    Table[i<->i+2, {i, 2, Length[rim]-4, 2 }],
    {Length[rim]-2 <-> 2},
    {Length[rim]-1 <-> 1}
  ]
, VertexCoordinates->rim]

rbonds = EdgeList[rbonds];
rpairs = {rbonds, Norm /@ getDir[rbonds, 2.0 rim]} // Transpose;

Начальные значения координат (для каждого шага по времени):

rim0 = 2.0 rim; rim1 = rim0; rim2 = rim0;
rimR = rim0;

mousePos   = {100.,0.};

Выведем их на экран, подпишемся на "анимацию", а также на движение мыши пользователя

EventHandler[
  showGeometry[rimR, rbonds, TransitionType->None, "Controls"->False, Epilog->{
    AnimationFrameListener[rimR // Offload, "Event"->"anim"], 
    Circle[mousePos // Offload, 0.5]
  }, PlotRange->{{-5,5}, {-5,5}}]
, {"mousemove" -> Function[xy,
    mousePos = xy;
]}]

Здесь мы подписываемся на цикл отрисовки браузера. Каждый раз когда браузер заканчивает рисовать кадр (обычно 1/60 сек), он отправляет событие "anim". В это время обработчик на стороне WL рассчитывает новое положение точек и отправляет его обратно браузеру по WebSocket. Важно отметить, что новое событие не будет инициировано до тех пор, пока не отработает старое. Таким образом обеспечивается максимально возможная частота кадров синхронизированная с обновлением экрана в текущей конфигурации системы (с учетом задержки сети, где находится сервер, клиент, производительности CPU, лагов JS и т.п.)

Теперь добавим также три кнопки управления анимацией, включая их обработчиков

Button["Start", 
 EventHandler["anim", With[{}, Do[With[{\[Delta]t = 0.006}, 
  rim0 = applyBond[rpairs, 1.0, 0.8][
    Clip[(2 rim1 - rim2) + Table[
       If[Max[Abs[i - mousePos]] < 2.0, {0,-1} - 7.0 (i - mousePos), {0,-1}]
    , {i, rim0}] Power[\[Delta]t,2], {-5,5}]
  ];

  rim2 = rim1;
  rim1 = rim0;
 ], {i, 10}];
  rimR = rim0;

 ]&]; 
 rimR = rim0;
]

Button["Stop", EventRemove["anim"]]

Button["Reset", With[{}, rim0 = 2.0 rim; rim1 = rim0; rim2 = rim0;
rimR = rim0;]]

Результат можете увидить ниже

Меш из изображения

Было бы здорово нарисовать что-то мышкой, и потом преобразовать в сетку вершин и связей автоматически

EventHandler[InputRaster[mask], (mask = #)&]

Нарисовали ножницы ✂️ хорошо. Теперь триангулируем и преобразуем в граф:

mesh = TriangulateMesh[
  ImageMesh[ImageResize[mask , 100] // ColorNegate]
, MaxCellMeasure->{"Area"->100}];

graph = MeshConnectivityGraph[mesh, 0] // IndexGraph
vertices = graph // GraphEmbedding;
edges = (List @@ #) &/@ (graph // EdgeList) ;

vertices = Map[Function[x, x - {0.4,-1.2}], 0.04 vertices];

Сделаем цикл анимации, как и раньше, и используем другие имена переменных

vertices0 =  vertices; vertices1 = vertices0; vertices2 = vertices0;
verticesR = vertices0;

vpairs = {edges, Norm /@ getDir[edges, vertices0]} // Transpose;

mousePos   = {100.,0.};

EventHandler[
  showGeometry[verticesR, edges, TransitionType->None, "Controls"->False, Epilog->{
    AnimationFrameListener[verticesR // Offload, "Event"->"anim2"], 
    Circle[mousePos // Offload, 0.5]
  }, PlotRange->{{-5,5}, {-5,5}}]
, {"mousemove" -> Function[xy,
    mousePos = xy;
]}]

Button["Start", 
 EventHandler["anim2", With[{}, Do[With[{\[Delta]t = 0.003}, 
  vertices0 = applyBond[vpairs, 1.0, 0.8][
    Clip[(2 vertices1 - vertices2) + Table[
       If[Max[Abs[i - mousePos]] < 2.0, {0,-1} - 7.0 (i - mousePos), {0,-1}]
    , {i, vertices0}] Power[\[Delta]t,2], {-5,5}]
  ];

  vertices2 = vertices1;
  vertices1 = vertices0;
 ], {i, 3}];
  verticesR = vertices0;

 ]&]; 
 verticesR = vertices0;
]

Button["Stop", EventRemove["anim2"]]

Button["Reset", With[{}, vertices0 =  vertices; vertices1 = vertices0; vertices2 = vertices0;
verticesR = vertices0;]]

Ура! У нас получилось ✨

Очевидно скриптовый язык плох в таких численных расчётах, и стоит поместить бОльшую часть логики в выражение Compile. Тогда вся тяжесть будет автоматически перенесена на сгенерированный и динамически подгруженный C-код. Но об этом в другой раз 🧙🏼‍♂️

Реализация базового метода Стёрмера-Верле

Набор точек

Процедурная генерация мешей

Обод

Меш из изображения

Публикации

Ближайшие события