ばねのシミュレーション

自力で考える

サスペンションのような硬いばねではなくヨーヨーをイメージする。そのヨーヨーは、手の位置に戻ってこようとするが、行きすぎてしまってなかなか手の位置に収まらない。

まずヨーヨーの本体と手の位置を定義する。

yoyo = 1.0
hand = 5.0

ヨーヨーの気持ちになると手の位置に戻りたいのだから 目的地の手 - 自分 で距離がわかる。

hand - yoyo  # => 4.0

それを自分に足せば、

yoyo += hand - yoyo  # => 5.0

となって手に移動できるので考え方としては間違っていない。間違っているとすれば瞬間移動してしまったところ。そこで距離を少しずつ足していく。

yoyo = 1.0
hand = 5.0
velocity = 0.0
a = 70.times.collect {
  yoyo += (hand - yoyo) * 0.2
  yoyo.round
}.join
a  # => "2233444445555555555555555555555555555555555555555555555555555555555555"

これで徐々に近づくようになった。しかし、ばねだったらぼよんぼよんしてほしい。そこで速度に足すように変更する。

yoyo = 1.0
hand = 5.0
velocity = 0.0
a = 70.times.collect {
  velocity += (hand - yoyo) * 0.2
  yoyo += velocity
  yoyo.round
}.join
a  # => "2357899875321123578998753211235789986532112457899865321124679998643111"

これで手の位置を中心にぼよんぼよんするようになった。これで完成でもいいけど、できれば収束させたい。そこで速度に減退率^[1]を掛けてみる。

yoyo = 1.0
hand = 5.0
velocity = 0.0
a = 70.times.collect {
  velocity += (hand - yoyo) * 0.2
  velocity *= 0.93
  yoyo += velocity
  yoyo.round
}.join
a  # => "2356788875433334556777655444445566666555444455566665555445555555555555"

これで手の位置に向けて収束するようになった。このように、専門的なことは置いといて、適当にやってもそれらしい動きになることがわかった。

リファクタリング

ここで、ばね業界について調べてみるとマジックナンバーに名前がついていることがわかった。その他の変数名も含めてリファクタリングした結果、

anchor = 5.0                    # 支点の位置
stiffness = 0.2                 # ばねの強さ (一言で言えば速度。大きいと速く動く)
damping = 0.93                  # 減退率 0..1 (大きいとぼよんぼよんが継続する)
mass = 1.0                      # 質量。大きいほど物体は動きにくくなる

location = 1.0                  # おもりの位置
velocity = 0.0                  # 現在の速度

a = 70.times.collect {
  force = (anchor - location) * stiffness
  acceleration = force / mass
  velocity += acceleration
  velocity *= damping
  location += velocity
  location.round
}.join
a  # => "2356788875433334556777655444445566666555444455566665555445555555555555"

となった。名前重要！

どの部分がばねの方程式になっているのか？

フックの法則のばねの方程式は F = -kx と定義されている。コードの (anchor - location) が -x に相当する。ここを間違えて逆に (location - anchor) としてしまうと x になってしまうので注意する。そして残った stiffness が k に相当する。

静止長を考慮する

ばねが伸びても縮んでもいない状態の長さを静止長(rest length)^[2]という。静止長を考慮するには支点を目標にするのではなく支点の少し手前を目標にする。その「少し手前」の幅が静止長に該当する。そうすると元のコード

ここは共通

anchor = 5.0                    # 支点の位置
location = 1.0                  # おもりの位置
stiffness = 0.2                 # ばねの強さ

静止長 0 の場合

force = anchor - location  # => 4.0
force *= stiffness         # => 0.8

は、

静止長を考慮した場合

rest_length = 1.0                       # 静止長
gap = anchor - location                     # => 4.0
stretch = gap.abs - rest_length             # => 3.0
force = (gap.positive? ? 1 : -1) * stretch  # => 3.0
force *= stiffness                          # => 0.6000000000000001

となる。このあたりスカラーだと符号の向きを気にかけないといけないので単にベクトルになっていた方が扱いやすい。ベクトルの場合は、

anchor = V[5.0, 0.0]            # 支点の位置
location = V[1.0, 0.0]          # おもりの位置
stiffness = 0.2                 # ばねの強さ

gap = anchor - location                  # => (4.0, 0.0)
stretch = gap.length - rest_length       # => 3.0
force = gap.normalize_or_zero * stretch  # => (3.0, 0.0)
force *= stiffness                       # => (0.6000000000000001, 0.0)

となって簡潔に書ける。

ヨーヨーのように静止長が 0 になる場合は平方根を求める処理が無駄になってしまうので、最適化するなら静止長が 0 かどうかで分岐するのがよさそう。

計算式まとめ

anchor = V[5.0, 0.0]            # 支点の位置

ばねの属性

stiffness = 0.01                # ばねの強さ (一言で言えば速度。大きいと速く動く)
rest_length = 1.0               # 静止長 (伸びても縮んでもない状態の長さ)

物体の属性

damping = 0.99                  # 減退率 0..1 (大きいとぼよんぼよんが継続する)
mass = 1.0                      # 質量 (大きいほど物体は動きにくくなる)
location = V[1.0, 0.0]          # 物体の位置
velocity = V.zero               # 現在の速度

物体の加速度を求める

gap = anchor - location                    # => (4.0, 0.0)
if rest_length == 0
  force = gap
else
  stretch = gap.length - rest_length       # => 3.0
  force = gap.normalize_or_zero * stretch  # => (3.0, 0.0)
end
force *= stiffness
acceleration = force / mass                # => (0.6000000000000001, 0.0)

加速度→速度→位置への反映

velocity += acceleration  # => (0.6000000000000001, 0.0)
velocity *= damping       # => (0.5580000000000002, 0.0)
location += velocity      # => (1.5580000000000003, 0.0)

最小限のコード

class App < Base
  def initialize
    super

    @stiffness = 0.02
    @damping   = 0.98
    @location  = window_wh * V.rand
    @velocity  = V.zero
    @mass      = 1.0
  end

  def update
    super

    force = (anchor - @location) * @stiffness
    acceleration = force / @mass
    @velocity = (@velocity + acceleration) * @damping
    @location += @velocity
  end

  def draw
    super

    point_draw(anchor)
    line_draw(anchor, @location)
    point_draw(@location)
    circle_draw(@location, 30)

    vputs "stiffness: #{@stiffness}"
    vputs "damping: #{@damping}"
  end

  def anchor
    mouse_v
  end

  show
end

パラメータの意味と考え方

いくつかあるパラメータのなかでもっとも影響があるのが「ばねの強さ」と「減退率」で、この調整が難しい。難しい原因は用語にある。「ばねの強さ」「減退率」と言われてもピンとこないからだ。そこで「速度」と「ぼよんぼよんする時間の長さ」と置き換えるとわかりやすい。単に「速度」と「時間」でもよい。

stiffness = 0.01                # ばねの強さ → 速度
damping   = 0.99                # 減退率 → ぼよんぼよんする時間の長さ → 時間

また、パラメータの初期値を適当にしてしまうとコードの不具合なのか区別がつかなくなるので初期値は上のようにしたい。

また、なんとなくだが、両方を足して 1.0 になるようにすると、あまり変な動きにはならないことに気づいた。たとえば速度を 0.02 にしたときは、時間を 0.98 にする。なので、調整しすぎてよくわからなくなってきたときには、いったん足して 1.0 を意識する。

密結合を改善する

• ばねの名称は「支点」「ばね」「おもり」に分けられる
• それをそのままオブジェクトと考えて問題ないだろうか？
• 少なくとも「おもり」が「ばね」に依存するのはおかしい
  • The Nature of Code では Spring と Bob (おもり) で分離されている
• 一方「支点」は「ばね」の一部だろうか？
• 「支点」は巨大な岩と考えることができる
• であれば「支点」も「ばね」に依存するのはおかしい

したがって「支点」「ばね」「おもり」は別々のオブジェクトになるはず。

そう考えて実装してみたところ、ごちゃごちゃになっていたパラメータが、どのオブジェクトの属性なのか明確になった。たとえば、stiffness や rest_length はばねの属性だが damping はおもりの属性になる。

ただ、数珠繋ぎにする場合など、オブジェクト自身がどこに繋っているかを把握していないと、管理するのが大変になることもわかった。

モデル化したコード

class Spring
  def initialize
    @stiffness = 0.02
    @rest_length = 80.0
  end

  def connect(a, b)
    from_to(a, b)
    from_to(b, a)
  end

  def from_to(from, to)
    if from.respond_to?(:acceleration_add)
      gap = to.location - from.location
      stretch = gap.length - @rest_length
      force = (gap.normalize_or_zero * stretch) * @stiffness
      from.acceleration_add(force)
    end
  end
end

class Anchor
  def location
    $app.mouse_v
  end

  def draw
    $app.point_draw(location)
  end
end

class Bob
  attr_accessor :location

  def initialize
    @damping = 0.98
    @location = V.rand * $app.window_wh
    @velocity = V.zero
    @mass = 1.0
    @acceleration = V.zero
  end

  def acceleration_add(force)
    @acceleration += force / @mass
  end

  def update
    @velocity += @acceleration
    @velocity *= @damping
    @location += @velocity
    @acceleration *= V.zero
  end

  def draw
    $app.point_draw(@location)
    $app.circle_draw(@location, 30)
  end
end

class App < Base
  def initialize
    super

    @anchor = Anchor.new

    @spring1 = Spring.new
    @bob1 = Bob.new

    @spring2 = Spring.new
    @bob2 = Bob.new
  end

  def update
    super

    @spring1.connect(@bob1, @anchor)
    @bob1.update

    @spring2.connect(@bob2, @bob1)
    @bob2.update
  end

  def draw
    super

    @bob1.draw
    line_draw(@anchor.location, @bob1.location)

    @bob2.draw
    line_draw(@bob1.location, @bob2.location)

    @anchor.draw
  end

  show
end

連結情報を保持するオブジェクトを用意してみる

ばねの「ばね」部分をラップした Connector クラスを新しく用意し、ばねの両端に接続するオブジェクトを保持するようにしてみた。が、あまりよくなったようには思えない。Nature of Code でもオブジェクト指向の方針に正解はないと書かれていて難しさを感じる。

Connector クラスで連結情報をまとめる

class Connector
  attr_accessor :head
  attr_accessor :tail

  def initialize
    @spring = Spring.new
    @head = []
    @tail = []
  end

  def update
    @head.product(@tail) do |a, b|
      @spring.connect(a, b)
    end
  end
end

class Spring
  def initialize
    @stiffness = 0.02
    @rest_length = 80.0
  end

  def connect(a, b)
    from_to(a, b)
    from_to(b, a)
  end

  def from_to(from, to)
    if from.respond_to?(:acceleration_add)
      gap = to.location - from.location
      stretch = gap.length - @rest_length
      force = (gap.normalize_or_zero * stretch) * @stiffness
      from.acceleration_add(force)
    end
  end
end

class Anchor
  def location
    $app.mouse_v
  end

  def draw
    $app.point_draw(location)
  end
end

class Bob
  attr_accessor :location

  def initialize
    @damping = 0.98
    @location = V.rand * $app.window_wh
    @velocity = V.zero
    @mass = 1.0
    @acceleration = V.zero
  end

  def acceleration_add(force)
    @acceleration += force / @mass
  end

  def update
    @velocity += @acceleration
    @velocity *= @damping
    @location += @velocity
    @acceleration *= V.zero
  end

  def draw
    $app.point_draw(@location)
    $app.circle_draw(@location, 30)
  end
end

class App < Base
  def initialize
    super

    @anchor = Anchor.new

    @spring1 = Spring.new
    @bob1 = Bob.new

    @spring2 = Spring.new
    @bob2 = Bob.new

    @connector1 = Connector.new
    @connector1.head << @anchor
    @connector1.tail << @bob1

    @connector2 = Connector.new
    @connector2.head << @bob1
    @connector2.tail << @bob2
  end

  def update
    super

    @connector1.update
    @connector2.update

    @bob1.update
    @bob2.update
  end

  def draw
    super

    @bob1.draw
    line_draw(@anchor.location, @bob1.location)

    @bob2.draw
    line_draw(@bob1.location, @bob2.location)

    @anchor.draw
  end

  show
end

バリエーション

計算方法は同じでも

• パラメータ調整
• ばねの連結
• 支点の有無

などによっていろいろな表現ができる。

ばねというより慣性がわずかに残る移動

オフストリング系ヨーヨー

ばねは両端のおもりの影響を受ける

ソロハム(失敗)

謎の生物

芋虫

夜店で売ってたばね

釣ったウナギ

実験用のコード

class Bob
  attr_accessor :location
  attr_accessor :velocity

  attr_accessor :parent
  attr_accessor :child

  delegate *[
    :gravity,
    :stiffness,
    :damping,
    :rest_length,
    :mass,
    :natural_mode,
    :gravity_mode,
    :mass_mode,
    :chain_mode,
    :anchor_mode,
    :wind_mode,
  ], to: "$app"

  def initialize(index)
    @index = index
    @location = $app.window_wh * V.rand
    @velocity = V.zero
    @acceleration = V.zero
    @radius = 30
  end

  def spring_update
    if anchor_mode && @index == 0
      spring_apply($app.mouse_v)
    end
    if chain_mode == 0
      # TODO: 重力がかかると子になるにつれて親との距離が半分になっていくるのは自然現象としてあっているのだろうか？
      spring_object_apply(parent)
      spring_object_apply(child)
    end
    if chain_mode == 1
      spring_object_apply(parent)
    end
    if chain_mode == 2
      spring_object_apply(child)
    end
  end

  def acceleration_add(force)
    if mass_mode
      @acceleration += force / mass
    else
      @acceleration += force
    end
  end

  def spring_object_apply(object)
    if object
      spring_apply(object.location)
    end
  end

  def spring_apply(anchor)
    gap = anchor - @location
    if natural_mode
      stretch = gap.length - rest_length
      force = gap.normalize_or_zero * stretch
    else
      force = gap
    end
    force *= stiffness
    acceleration_add(force)
  end

  def update
    spring_update

    if wind_mode
      acceleration_add(V.x * 1.0)
    end

    @velocity += @acceleration
    @velocity *= damping

    if gravity_mode
      @velocity += gravity
    end

    @location.replace(@location + @velocity)
    @acceleration = V.zero
  end

  def draw
    if child
      $app.line_draw(@location, child.location)
    end
    $app.point_draw(@location)
    $app.circle_draw(@location, @radius)
    if anchor_mode && @index == 0
      $app.line_draw($app.mouse_v, @location)
    end
    if false
      if parent
        $app.vputs "bob: [#{@index}] #{@location.round(2)} #{(parent.location - location).length}"
      end
    end
  end
end

class ParticleSystem < Array
  def update
    each(&:update)
  end

  def draw
    each(&:draw)
  end
end

class App < Base
  attr_accessor :gravity
  attr_accessor :stiffness
  attr_accessor :damping
  attr_accessor :rest_length
  attr_accessor :mass

  attr_accessor :natural_mode
  attr_accessor :gravity_mode
  attr_accessor :mass_mode
  attr_accessor :chain_mode
  attr_accessor :anchor_mode
  attr_accessor :wind_mode
  attr_accessor :muscle_mode
  attr_accessor :ring_mode

  def initialize
    super

    @wind_mode = false
    preset1
  end

  def ps_init
    @ps = ParticleSystem.new
    ps_add(@bob_count)
  end

  def ps_add(n = 1)
    n.times do
      a = @ps.last
      b = Bob.new(@ps.count)
      if a
        a.child = b
        b.parent = a
      end
      @ps.concat([b])
    end
    @points = @ps.collect(&:location)
    ring_mode_update
  end

  def ps_remove(n = 1)
    n.times { @ps.pop }
    @points = @ps.collect(&:location)
    ring_mode_update
  end

  def ring_mode_update
    unless @ps.empty?
      a = @ps.first
      b = @ps.last
      if @ring_mode
        a.parent = b
        b.child = a
      else
        a.parent = nil
        b.child = nil
      end
    end
  end

  def ring_mode_toggle
    @ring_mode ^= true
    ring_mode_update
  end

  # 基本
  def preset1
    @natural_mode = false
    @gravity_mode = false
    @anchor_mode = true
    @ring_mode = false
    @muscle_mode = false

    @mass_mode = false
    @chain_mode = 0

    @bob_count = 1
    ps_init

    @gravity = V[0.0, 1.0]
    @stiffness = 0.02
    @damping = 0.96
    @mass = 3.0
    @rest_length = 100.0
  end

  # 速い移動
  def preset2
    preset1

    @stiffness = 0.2
    @damping  = 0.7
  end

  # ヨーヨー 4A
  def preset3
    preset1

    @gravity_mode = true
    @anchor_mode = true
    @stiffness = 0.01
    @damping  = 0.98
    @gravity  = V[0.0, 0.7]
  end

  # ヨーヨー 5A
  def preset4
    preset1

    @anchor_mode = false
    @natural_mode = true
    @ring_mode = false
    @muscle_mode = false

    @bob_count = 2
    ps_init
  end

  # ソロハム
  def preset5
    preset1

    @anchor_mode = true
    @natural_mode = true
    @ring_mode = true
    @muscle_mode = false
    @gravity_mode = true

    @stiffness = 0.02
    @damping = 0.98
    @rest_length = 100.0
    @gravity = V[0.0, 0.5]
    @bob_count = 3
    ps_init
  end

  # 謎の生物
  def preset6
    preset1

    @anchor_mode = false
    @natural_mode = true
    @ring_mode = true
    @muscle_mode = true
    @gravity_mode = false

    @stiffness = 0.02
    @damping = 0.96
    @rest_length = 70.0
    @bob_count = 3
    ps_init
  end

  # 芋虫
  def preset7
    preset1

    @natural_mode = true
    @gravity_mode = false
    @anchor_mode = true
    @ring_mode = false
    @muscle_mode = false

    @stiffness = 0.6
    @damping  = 0.6
    @rest_length = 25.0
    @bob_count = 20
    ps_init
  end

  # 夜店で売ってたばね
  def preset8
    preset1

    @natural_mode = true
    @gravity_mode = false
    @anchor_mode = true
    @ring_mode = false
    @muscle_mode = false

    @stiffness = 0.7
    @damping  = 0.6
    @rest_length = 5.0
    @bob_count = 30
    ps_init
  end

  # 釣ったウナギ
  def preset9
    @natural_mode = true
    @gravity_mode = true
    @anchor_mode = true
    @ring_mode = false
    @muscle_mode = false

    @stiffness = 0.45
    @damping = 0.7
    @rest_length = 3.0
    @gravity = V[0.0, 0.6]
    @bob_count = 20
    ps_init
  end

  def button_down(id)
    super

    case id
    when Gosu::KB_1
      preset1
    when Gosu::KB_2
      preset2
    when Gosu::KB_3
      preset3
    when Gosu::KB_4
      preset4
    when Gosu::KB_5
      preset5
    when Gosu::KB_6
      preset6
    when Gosu::KB_7
      preset7
    when Gosu::KB_8
      preset8
    when Gosu::KB_9
      preset9
    when Gosu::KB_N
      @natural_mode ^= true
    when Gosu::KB_G
      @gravity_mode ^= true
    when Gosu::KB_S
      @mass_mode ^= true
    when Gosu::KB_W
      @wind_mode ^= true
    when Gosu::KB_M
      @muscle_mode ^= true
    when Gosu::KB_C
      @chain_mode = @chain_mode.next.modulo(3)
    when Gosu::KB_R
      ring_mode_toggle
    when Gosu::KB_A
      @anchor_mode ^= true
    when Gosu::KB_LEFT, Gosu::KB_RIGHT
      @stiffness += (id == Gosu::KB_LEFT ? -1 : 1) * 0.01
    when Gosu::KB_UP, Gosu::KB_DOWN
      @damping += (id == Gosu::KB_UP ? 1 : -1) * 0.01
      @damping = @damping.clamp(0.0..1.0)
    when Gosu::KB_J, Gosu::KB_K
      @rest_length += (id == Gosu::KB_J ? 1 : -1) * 1.0
    when Gosu::KB_H, Gosu::KB_L
      @gravity += (id == Gosu::KB_H ? -1 : 1) * 0.1 * V.y
    when Gosu::KB_RETURN
      @bob_count += 1
      ps_add
    when Gosu::KB_BACKSPACE
      @bob_count += -1
      ps_remove
    end
  end

  def draw
    super

    vputs "[←][→] ばね係数(stiffness): #{@stiffness.round(3)}"
    vputs "[↑][↓] 減退度(damping): #{@damping.round(2)}"
    vputs "[j][k] rest_length: #{@rest_length.round(2)}"
    vputs "[h][l] gravity: #{@gravity.round(2)}"
    vputs "[RET][BS] bob_count: #{@bob_count}"
    vputs
    vputs "[g] gravity #{@gravity_mode}"
    vputs "[n] rest_length #{@natural_mode}"
    vputs "[s] mass #{@mass} #{@mass_mode}"
    vputs "[m] muscle: #{@muscle_mode}"
    vputs "[a] anchor_mode: #{@anchor_mode}"
    vputs "[r] ring_mode: #{@ring_mode}"
    vputs "[w] wind_mode: #{@wind_mode}"
    vputs "[c] chain_mode: #{@chain_mode} (experimental)"
    vputs
    vputs "Preset"
    vputs "[1] 基本"
    vputs "[2] 速い移動"
    vputs "[3] ヨーヨー 1A / 4A (aで離す)"
    vputs "[4] ヨーヨー 5A"
    vputs "[5] ソロハム"
    vputs "[6] 謎の生物"
    vputs "[7] 芋虫"
    vputs "[8] 夜店で売ってたばね"
    vputs "[9] 釣ったウナギ"
    vputs

    @ps.update
    @ps.draw
    point_draw(mouse_v)
  end

  def rest_length
    if @muscle_mode
      @rest_length * 0.5 + Math.sin((@frame_count / (60.0 * 1.0)).turn_to_rad) * @rest_length * 0.5
    else
      @rest_length
    end
  end

  def window_size_default
    V[1200, 800]
  end

  def needs_cursor?
    false
  end

  show
end

参照

• https://natureofcode.com/book/chapter-2-forces/
• https://yoppa.org/proga10/1271.html
• https://ja.wikipedia.org/wiki/フックの法則

脚注

1. 正確には減退しない率を掛けているのだから非減退率というべき？

2. 自然長(natural length)という呼び方もあるようだが、ばねに限定せず広い範囲で使われる名前のようなので、静止長(rest length)の方を使うことにする