Реализация полуявного обратного Эйлера в системе масс-пружин 1-DOF - PullRequest
Новая
       12

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

6 голосов
/ 09 октября 2010

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

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

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

Подводя итог, это приводит к

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

Применяя, например, Явный Эйлер, можно получить следующее:

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

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

1 Ответ

14 голосов
/ 09 октября 2010

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

Неявный Эйлер будет иметь (давайте назовем эти eqn (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 '( (x, v) 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

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

...