Реализация полунеявного обратного Эйлера в системе масса-пружина с 1 степенью свободы

У меня есть простая (массовая) пружинная система с двумя точками, соединенными пружиной. Одна точка закреплена на потолке, поэтому я хочу рассчитать положение второй точки численным методом. Итак, в основном я получаю положение второй точки и ее скорость и хочу знать, как эти два значения обновляются после одного временного шага.

На точку действуют следующие силы:

  • Гравитационная сила, заданная -g * m
  • Сила пружины, определяемая как k * (l - L), где k — жесткость, l — текущая длина, а L — начальная длина.
  • Демпфирующая сила, заданная -d * v

В сумме это приводит к

  • F = -g * m + k * (l - L)
  • Fd = -d * v

Применяя, например, Explicit Euler, можно вывести следующее:

  • newPos = oldPos + dt * oldVelocity
  • newVelocity = oldVelocity + dt * (F + Fd) / m, используя F = m * a.

Однако теперь я хочу использовать полунеявный обратный Эйлер, но не могу точно понять, откуда брать якобианы и т. д.


interpolation algorithm simulation
person Etan    schedule 09.10.2010    source источник
comment
удален тег spring, так как он связан с фреймворком java spring   -  person Sean Patrick Floyd    schedule 11.10.2010

Ответы (1)


arrow_upward
14
arrow_downward

Поэтому, вероятно, проще всего увидеть, как это происходит, если сначала рассмотреть полностью неявный метод, а затем перейти к полунеявному.

Неявный Эйлер имел бы (назовем это уравнением (1)):

newPos = oldPos + dt * newVelocity
newVelocity = oldVelocity + dt * (-g * m + k*(newPos - L) - d*newVelocity)/m

А пока давайте просто измерим позиции относительно L, чтобы избавиться от термина -kL. Переставляя, мы получаем

(newPos, newVelocity) - dt * (newVelocity, k/m newPos - d/m newVelocity) = (oldPos, oldVelocity - g*dt)

и представить это в матричной форме

((1,-dt),(k/m, 1 - d/m)).(newPos, newVelocity) = (oldPos, oldVelocity -g*dt)

\left ( \begin{array}{cc} 1 & -\Delta t \\frac{k}{m} & 1 - \frac{d}{m}\end{array} \right ) \left ( \begin{array}{c} x^\mathrm{new} \ v^\mathrm{new} \end{array} \right ) = \left ( \begin{array}{c} x^\mathrm{old} \ v^\mathrm{old} - g \Delta t\end{array} \right )

Где вы знаете все в матрице, и все в RHS, и вам просто нужно решить для вектора (newPos, newVelocity). Вы можете сделать это с помощью любого решателя Ax=b (исключение Гаусса вручную работает в этом простом случае). Но поскольку вы упоминаете якобианцев, вы, вероятно, хотите решить эту проблему с помощью итерации Ньютона-Рафсона или чего-то подобного.

В этом случае вы, по сути, хотите решить нули уравнения

((1,-dt),(k/m, 1-d/m)).(newPos, newVelocity) - (oldPos, oldVelocity -g*dt) = 0

то есть f(newPos, newVelocity) = (0,0). У вас есть предыдущее значение для использования в качестве начального предположения (oldPos, oldVelocity). Теперь вы просто хотите повторить

(x,v)n+1 = (x,v)n + f((x,v)n)/f'( (х,в)n)

пока вы не получите достаточно хороший ответ. Здесь,

f(newPos,newVel) = ((1,-dt),(k/m, 1-d/m)).(newPos, newVelocity) - (oldPos, oldVelocity -g*dt)

а f'(newPos, newVel) — якобиан, соответствующий матрице

((1,-dt),(k/m, 1-d/m))

Прохождение процесса для полунеявного является таким же, но немного проще - не все члены RHS в уравнениях (1) являются новыми величинами. Как это обычно делается, это

newPos = oldPos + dt * newVelocity
newVelocity = oldVelocity + dt * (-g * m + k*oldPos - d*newVelocity)/m

например, скорость зависит от старого временного значения положения, а положение — от нового временного значения скорости. (Это очень похоже на интеграцию «чехарда»..) Вы должны быть в состоянии довольно легко выполнить описанные выше шаги с этим немного другим набором уравнений. По сути, член k/m в приведенной выше матрице исчезает.

person Jonathan Dursi    schedule 09.10.2010
comment
Между прочим, Гилберт Стрэнг (который имеет большое значение) читал в Массачусетском технологическом институте курс именно по таким вещам, который есть на OpenCourseWare, и я думаю, что лекции просто превосходны: ocw.mit.edu/courses/mathematics/ - person Jonathan Dursi; 09.10.2010
comment
Вы можете использовать разметку TeX. Если я не ошибаюсь, вы начинаете и заканчиваете знаком доллара, чтобы обозначить TeX. - person JoshD; 10.10.2010
comment
Похоже, это только по математике и, возможно, латексным переполнениям стека =(. Но при поиске нашел mathurl.com и латексную лабораторию Google; классная штука. - person Jonathan Dursi; 10.10.2010