Физика вращения 3д тел

Думаю, что ошибка в update

К сожалению да, в примере ошибка, которую в коде я давно исправил, а там забыл. В коде сейчас вроде бы корректно.

На момент построения графиков было так:

class State3d(val transform: Transform3d, val velocity: Velocity3d):
  def updated(s: State3dDerivative, dt: Double): State3d =
    State3d(
      transform.copy().update(s.velocity, dt),
      velocity.copy().update(s.acceleration, dt)
    )

Я явно копирую transform и velocity, state3d не меняется, создаётся новая копия.

class Transform3d(val position: Vector3d, val rotation: Quaternion):
  def update(velocity: Velocity3d, dt: Double): Transform3d =
    position.madd(velocity.linear, dt)
    val dRot = Quaternion().setFromExp(velocity.angular * (0.5 * dt))
    dRot *> rotation
    rotation.normalize()
    this

class Velocity3d(val linear: Vector3d, val angular: Vector3d):
  def update(acceleration: Acceleration3d, dt: Double): Velocity3d =
    linear.madd(acceleration.linear, dt)
    angular.madd(acceleration.angular, dt)
    this

Я посмотрел Ваш пример, но в нём производятся вычисления в локальной системе отсчёта и вообще не затрагивается вопрос ориентации тела в пространстве. Попробую его воспроизвести и прям по шагам сравнить вычисления со своими, но не уверен, что быстро найду ошибку.

Большое спасибо за пример!

В Вашем примере в качестве ошибки берётся (|L| — |L0|) / |L0|. Если в моём коде считать ошибку так, то она тоже получется порядка 10^-8.
Но если считать ошибку как |L — L0|/|L0|, то получается, что L отклоняется в сторону и получается ошибка порядка 10^-4.

Я повторил Ваши рассчёты в локальной СО для скорости и ускорений. Они работают с такой же точностью: 10^8 при dt =0.01. И для E, и для модуля L.

После этого я попробовал добавить обновление ориентации. Рассчёты в локальной СО от неё никак не зависят, но станет можно повернуть L из локальной СО в глобальную. В глобальной СО L не должен отклоняться. К сожалению, он отклоняется, и итоговая ошибка тоже получается порядка 10^-4 (как и у меня). По крайней мере, у меня лучше не получилось :(

У double мантисса 52 бита, это примерно 16 значащих цифр. Я думаю, такой точности с запасом хватит почти для всего, кроме каких-то спецефических случаев.