Document

The Virtual Reality Society of Japan
インタラクティブ技術に必要な
数学と物理
東京工業大学
長谷川晶一
The Virtual Reality Society of Japan
自己紹介
長谷川晶一
VRをはじめて11年．学部２年からVRコンテストに参加
コンテストで忙しくて授業に出ない日々
教養の数学・物理は，ぎりぎりで単位取得
VRコンテスト→東工大佐藤・小池研
（力覚インタフェースSPIDARの研究室）
最近は動力学シミュレーションも研究
The Virtual Reality Society of Japan
VRシステムの中の数学と物理
VRシステム
提示装置
計測装置
大画面
HMD
ポヒマス
加速度センサ
コンピュータ
3DCGレンダリング
シミュレーション
The Virtual Reality Society of Japan
今日の目標
数学がVR役立つということを納得してもらいたい．
VRシステムを作るために必要な最小限の数学と
物理の話をする．
あとで，自分の問題を解くために，教科書のどこを見れば
よいか分かるように．
今日の話は(あたりまえですが)不完全です．いつか，
それに気づいてほしいです．
The Virtual Reality Society of Japan
話の進行
トピック
アプリケーション
位置の数値化
位置・角度センサ
座標系の変換と行列
逆変換
回転変換
行列のランク
3DCG
ロボット
積分・微分と
シミュレーション
シミュレーション
物理シミュレーション
運動方程式（質点・剛体）
安定性と行列のn乗制御
教科書
線形代数
微分積分
力学
線形代数
The Virtual Reality Society of Japan
位置と数値
ポヒマスセンサという位置センサ
トランスミッタからみた，レシーバ位置を計測
トランスミッタ
レシーバ
The Virtual Reality Society of Japan
位置と数値
数値で表さなければならない
コンピュータ
トランスミッタ
レシーバ
The Virtual Reality Society of Japan
位置と数値
P
O


p
レシーバの位置(P)は，トランスミッタ(O)から
矢印pだけ進んだところ
矢印pを数値で表したい．
The Virtual Reality Society of Japan
位置と数値
P
O
p
ey
ez



レシーバの位置(P)は，トランスミッタ(O)から
矢印pだけ進んだところ
矢印pを数値で表したい．
向きと長さの基準が必要：基底(基準の矢印)
矢印p は，ex ３つ分+ ey ２つ分+ ez １つ分
これを，p=3ex + 2ey + 1ez と書く
ex
The Virtual Reality Society of Japan
位置と数値

レシーバの位置(P)は，トランスミッタ(O)から
矢印p =3ex + 2ey + 1ezだけ進んだところ
P
O
ez

ey
p
ex
原点O，基底(ex ey ez ) が決まると，
3
点Pは3つの数値 2 で表せる．
1

(原点O，基底(ex ey ez ) )のことを座標系といいます．
The Virtual Reality Society of Japan
位置と数値
ey
p=3ex + 2ey + 1ez
基底のおかげで
矢印が数値で表せるようになった
ex
ベクトル
3
p= 2
1
ez
基底は(ex ,ey ,ez )
基底を忘れてはならないのだけど．．．
面倒なので省略します．
ついでにいうと
1
0
基底ベクトル ex ey ez も
ex = 0 ey = 1
矢印だから数値で書けます．
0
0
それにも，もちろん，基底(ex ,ey ,ez )が必要です．
0
ez = 0
1
The Virtual Reality Society of Japan
座標系
手の位置の計測
トランスミッタを目の位置に置くと
ey
トランスミッタO
ez ph
ex
Ph
レシーバ
ph = 2ex -1ey +3ez
2
-1
3
目にくっついた座標系で，手の位置が数値になる．
The Virtual Reality Society of Japan
座標系の変換
机から見た手の位置
ポヒマスが2個あると
pe = 2wx +1wy +1wz
ey
wy
pe e
z
wx q
wz
机にくっついた座標系wで
目の位置が数値になる
ph
ex
ph = 2ex -1ey +3ez
目にくっついた座標系eで
手の位置が数値になる
2
1
1
2
-1
3
q = pe+ph ：机から見た手の位置
机にくっついた座標系wでの
手の位置を数値にすると？
The Virtual Reality Society of Japan
向きの計測
ポヒマスセンサは実は向きも測れる
wy
トランスミッタ
wx
wz
コンピュータ
ey
レシーバ
ez
ex
位置pと同様に向き(ex ey ez )も数値で返す
ey wy
0.5 0 -0.9
(ex ey ez )=( 0 1 0 )
0.9 0 0.5
基底は(wx wy wz )
ez
ex
wz
wx
The Virtual Reality Society of Japan
基底の変換（取替え）
ポヒマスの返す数値から
ex = 0.5 wx + 0 wy + 0.9 wz
ey = 0 wx + 1 wy + 0 wz
ez =-0.9 wx + 0 wy+ 0.5 wz
ey wy
wx
ez
だから，eをwで書き直せて，
wz
p
h
ph = 2ex -1ey +3ez
= 2(0.5 wx + 0.9 wz )-1wy +3(-0.9 wx +0.5 wz )
= (1-2.7)wx -1wy +(1.8+1.5)wz
= -1.7 wx -1wy +3.3 wz
q=pe+ ph =
2
-1.7
0.3
1 + -1 = 0
1
3.3
4.3
ex
机にくっついた座標系での
手の位置の数値
The Virtual Reality Society of Japan
基底の変換と行列
簡単に書きたい
x
と書くことにする．
y
ph = 2ex -1ey +3ez
ax+byを(a b)
2
ex = 0.5 wx + 0.9 wz
ph = (ex ey ez ) -1
ey = w y
ez =-0.9 wx + 0.5 wz なので，
3
0.5
0
-0.9
= (wx wy wz ) 0 (wx wy wz ) 1 (wx wy wz ) 0
0.9
0
0.5
2
-1
3
-1.7
= (wx wy wz ) -1
3.3
基底を数値で表したもの→行列
0.5 0 -0.9
= (wx wy wz ) 0 1 0
0.9 0 0.5
2
-1
3
The Virtual Reality Society of Japan
基底の変換と行列のまとめ
ey wy
phは，2ex +1ey +3ez
2
phは， 1
3
これをwで表したい
0.5
(ex ey ez )は，( 0
0.9
0
1
0
wx
ez
-0.9
0 )
0.5
なので，phは，
0.5 0 -0.9
ph= 0 1 0
0.9 0 0.5
行列
wz
ex
p 2
1
3
-1.7
2
1 = -1
3.3
3
行列とは基底を数値で書いて並べたもの
The Virtual Reality Society of Japan
いろいろな変換
2DCGの変形
wy
線の描画→ 始点と終点を数値で指定
基底：画面上のwx wy
基底を変えると変形する．
wx
ey
ex
ey
ey
(ex ey )=(
0 0
)
0 0.5
(ex ey )=( 2
0
(ex ey )=(
0 )
1
1 1
)
1 1
ey
ex
ex
(ex ey )=( 1 0 )
0 1
ex
(ex ey )=( 0
1
1 )
0
(ex ey )=(
1 -1
)
1 -1
(ex ey )=( 1 0 )
1 1
The Virtual Reality Society of Japan
回転変換
 回転行列
 cos
R( )  
 sin 
0
wy 1
ey
 sin  

cos 
というのがあったような
wy
1
ex 0


wx 1
0
ey
 cos 

e x  
sin



sin
  sin  

e y  
wx
cos



ex

cos
回転変換→回転した基底での数値をもとの基底での数値に変換する．
The Virtual Reality Society of Japan
逆変換
机の上の1台のトランスミッタで
目の座標系で表した手の位置を計測したい
ey
wy
O
wz
pe
q ez
wx
ph
ex
ポヒマスの出力：
pe q ex ey ezの数値
基底は(wx wy wz )
知りたいもの：
ph1
phの数値
基底は(ex ey ez ) ph2
ph3
ph1w
ph=q-pe なので，
基底(wx wy wz )でのph数値 ph2w は
ph3w
わかる
The Virtual Reality Society of Japan
逆変換
(wx wy wz)を基底にしたex ey ezの数値は
ポヒマスから得られるので，
ph1
ph1w
ph2w = ex ey ez ph2
ph3
ph3w
もし，逆に(ex ey ez) を基底としたwx wy wzの数値がわかれば，
ph1
ph1w
ph2 = wx wy wz ph2w
ph3
ph3w
ph1
なので，ph2 の数値がわかる．
ph3
The Virtual Reality Society of Japan
逆変換
ph1w
ph2w =
ph3w
ex ey ez
ph1
ph2
ph3
ph1w
ph2w =
ph3w
ex ey ez
wx wy wz
ph1w
ph2w =
ph3w
1 0 0
0 1 0
0 0 1
ph1w
ph2w
ph3w
ex ey ez
逆変換
逆行列
ph1
ph2 =
ph3
ph1w
ph2w
ph3w
wx wy wz
=
1 0 0
0 1 0
0 0 1
-1
ex ey ez =
wx wy wz
と書く
ph1w
wx wy wz ph2w
ph3w
The Virtual Reality Society of Japan
逆行列の求め方
ex ey ez
x2倍
z=z- √3/2 x
1/2
0
√3/2
x2倍 1
0
√3/2
1
0
z=z- √3/2 x 0
0 - √3/2
1 0
0 1/2
0 - √3
1 0
0 1/2
0 - √3
1 0
0 2
wx wy wz
1
0
0
2
0
0
2
0
- √3
0
1
0
0
1
0
0
1
0
0
0
1
0
0
1
0
0
1
The Virtual Reality Society of Japan
逆行列の求め方
z=z- √3/2 x
1
0
z=z- √3/2 x 0
0 - √3
1 0
0 2
2 0 0
0 1 0
- √3 0 1
1
0
0
0 - √3
1 0
0 1
2
0 0
0
1 0
- √3/2 0 1/2
1
0
0
0
1
0
0
0
1
1/2 0 √3/2
0 1 0
- √3/2 0 1/2
ex ey ez
wx wy wz
1/2倍
z 1/2倍
x=x+ √3z
x=x+ √3z
The Virtual Reality Society of Japan
3DCGでの変換の利用
ポリゴン・メッシュとは，3角形の頂点座標
1 0 1
W: ( 0 1 0 )とO
0 0 1
M:(mx my mz )とmp
ワールド座標系
モデリングするときの座標系
V:(vx vy vz )とvp 視点座標系（視点の位置と向き）
my
p3
p2
mx
wy
vy
vx
mp mz
O
wz
wx
視点
p1
vp
vz
The Virtual Reality Society of Japan
３DCGの変換
 モデリングしたときの頂点位置p1 p2 p3は，M座標系での数値になってい
る．数値：p1m p2m p3m
 これをV座標系での数値：p1ｖ p2ｖ p3ｖにできれば画面に表示できる．
my
v
p3
p2
wy
vy
mp m
z
O
wz
y
mx
p2
vx
視点
p1
vp
wx
p1w= mp + Mp1m
p1w=vp + Vp1v
p1v = V-1 (p1w -vp)
p3
vz
vz
vx
p１
p1v =V-1 (mp +Mp1m-vp )
p1v x p2v x
あとは， p ， p
…に線を引けばよい
2v y
1v y
The Virtual Reality Society of Japan
逆変換のない場合
戻せる
ey
ey
ey
ex
ex
ey
ex
1
(ex ey )=( 0
ex
0
1 )
(ex ey )=(
2
0
0
)
1
ey
(ex ey )=(
0
1
1
)
0
1
(ex ey )=( 1
戻せない
(ex ey )=(
0 0
)
0 0.5
(ex ey )=(
ex
1 1
)
1 1
つぶれてしまうと元に戻せない
(ex ey )=(
1 -1
)
1 -1
0
1 )
The Virtual Reality Society of Japan
行列のランク
ロボットの関節の速度と手先の速度
e1: 1が1のときの手先の速度
v
  r1 sin 1  r2 sin  2 


r
cos


r
cos

1
2
2 
 1
e1
wy
r2
2  2
e2
wx O
e2: 2が1のときの手先の速度
1 1
r1
  r2 sin  2 


 r2 cos  2 
ロボットなので，
r1, r2, 1, 2,はわかる．
1 , 2が設定できる．
手先に速度vを出したい
The Virtual Reality Society of Japan
行列のランク
1 , 2がわかっているとき，手先の速度は
v  1e1  2e 2
 1 
 e1 e 2  
 2 
  r1 sin 1  r2 sin  2
 
 r1 cos1  r2 cos 2
 r2 sin  2   1 
  
r2 cos 2   2 
となる．
vができたのだから，逆変換 v  をすれば，
vを実現するがわかる
 1    r1 sin 1  r2 sin  2
   
 2   r1 cos1  r2 cos 2
 r2 sin  2  1
 v
r2 cos 2 
逆変換できない場合は？
The Virtual Reality Society of Japan
行列のランク
e1
r1=1, r2=1, 1=30°, 2 =90°のとき
e1
  r1 sin 1  r2 sin  2

e 2  
 r1 cos1  r2 cos 2
 r2 sin  2    1.5  1
  

r2 cos 2   3 / 2 0 
基底が2つある．
e2
逆変換あり．どんなvでも出せる．
r1=1, r2=0, 1=90°, 2 =90°のとき
e1
 1 0 

e 2   
 0 0
基底が１つしかない．
e2
逆変換なし．vは直線上だけ
r1=1, r2=1, 1=30°, 2 =30°のとき
 1 1/ 2 

e1 e2   
3 / 2
 3
実は基底が１つしかない．
e1
e1
e2
逆変換なし．vは直線上だけ
The Virtual Reality Society of Japan
行列のランク
r1=0, r2=0, 1=30°, 2 =90°のとき
e1
基底がない 0個．
 0 0

e2   
 0 0
e1 e 2
逆変換なし．vはまったく出せない
基底の数のことを行列のランクといいます．
・平行な基底は数えません．
・長さが0の基底は数えません．
2次元なら2つ，3次元なら3つ，4次元なら4つ ….
基底がないと，つぶれてしまい，逆変換できません
The Virtual Reality Society of Japan
冗長マニピュレータと長方形行列
 v x   e1x
  
v   e
 y   1y
e2 x
e2 y
 1 
e3 x   
  2 
e3 y   
 3 
e1
e2
wy
基底は3つあるように見えるが，
e1
e2 = 0.5e1 +0.8e3
のように，ほかの2つで表せる
e3
wx
速度vは2次元なのに，は3つ
3 3
O
1 1
2  2
冗長マニピュレータ
0.5
e2
e3
0.8
関節2を動かさなくても，関節1と3を動かせば，
どんな速度vも出せる．
逆行列はない．答えが無数にある．
The Virtual Reality Society of Japan
シミュレーション
 シミュレーション
 現実のものの動きを計算機の中に再現すること
 そのために，現実世界の法則を見つけ，それを再現する．
 例：車の動きのシミュレーション
 「アクセルを踏むと加速する」→法則
アクセル:原因
時刻t=0で，
0m/sだった
シミュレーション
速度
？
速度:結果
アクセルの踏み具合
時間t
The Virtual Reality Society of Japan
シミュレーション
アクセル:原因
車の動き:結果
シミュレーション
「アクセルを踏むと車は加速する」法則を式で書くと
ちょっと後（t+t）の速度= 今（t）の速度+ 今のアクセル×踏んだ時間
v(t  t )  v(t )  a(t )t
今のアクセル
t t +t
tが小さければ a(t) と a(t+t)はだいたい等しい
tが小さければ小さいほど誤差が少ない．
The Virtual Reality Society of Japan
数と関数
a(t )は， atと似ている．
ただ at は a0とか a1 のように，右下は整数だった．
数がならんでいると思
っていた．
a1
a(t )は， a(0)とか a(1)だけでなく a(0.1)とか a(1 / 2)とか
a( 3 )とかも考える．
a2 a3 a4
a(t)
a0
0
ここまでは，aといったら変数だったので値は1つだった．
a(t)は，a0 a１a2 …を含んでいる．
だけど，これからは，a(t)と書かないでaと書く．こともある．
この場合，aは関数．たくさんの値を持っていて，tごとに
値a(t)がある．
t
The Virtual Reality Society of Japan
シミュレーション
法則 v(t  t )  v(t )  a(t )t
分かっていること
最初の状態 v(0)  0m / s
v(0.1)  v(0)  a (0)0.1  0m / s
v(0.2)  v(0.1)  a (0.1)0.1  0.005m / s
0.05
その後の v(0.3)  v(0.2)  a0.1(0.2)0.1  0.015m / s
v(0.4)  v(0.3)  a (0.3)0.1  0.03m / s
経過
0.15
v(0.5)  v(0.4)  a (0.4)0.1  0.05m / s
1
2
0.2
．
．
v(0.9)  v(0.8) ．a(0.8)0.1  0.18m / s
0.4
v(1.0)  v(0.9)  a(0.9)0.1  0.225m / s
0.45
これをプログラムでやるのが，コンピュータシミュレーション
The Virtual Reality Society of Japan
シミュレーションと積分
v(t  t )  v(t )  a(t )t
v(t)  v(t  t )  a(t )t
 v(t  2t )  a(t  t )t  a(t )t
 v(t  3t )  a(t  2t )t  a(t  t )t  a(t )t

 v(0)  (a(t )  a(2t )  a(3t )    a(t ))t
t / t
v(t )  v(0)   a (it )t
i 0
誤差をなくすためtを0に近づける
t
積分
 v(0)   a( )d
 0
The Virtual Reality Society of Japan
シミュレーションと積分
シミュレーションではtずつ時間を進めて
ほしい時刻まで繰り返し計算する．
v(t  t )  v(t )  a(t )t
積分はその結果のこと
t / t
v(t )  v(0)   a(it )t
i 0
誤差をなくすためtを0に近づける
t
 v(0)   a( )d
 0
The Virtual Reality Society of Japan
シミュレーションと微分
車の速度v(t)の記録がある．
アクセルの踏み具合a(t)を知りたい
v (t  t )  v (t )  a (t ) t
v (t  t )  v (t )
a (t ) 
t
v(t)
0.275
0.225
v(1  0.1)  v(1)
a(1) 
 0.5m / s 2
0.1s
1.0 1.1
v(t)がわかっているので，
ちょっと後の値と今の値を引き算すれば出てくる．
The Virtual Reality Society of Japan
シミュレーションと微分
v (t  t )  v (t )
a (t ) 
t
誤差をなくすため，tを0に近づける
dv (t )

微分
dt
dv (t ) vをtで微分する．
vのt微分
dt
dv (t ) 微分が入った方程式
a (t ) 
→微分方程式と呼ぶ
dt
The Virtual Reality Society of Japan
微分と積分
v(t  t )  v(t )  a(t )t
(ただしtは十分小さい)
これが元の式. a(t)とv(t)の関係を表している．
v(t  t )  v(t ) dv (t )
a(t ) 

t
dt
t / t
微分
v(t )  v(0)   a(it )t  v(0)  
i 0
位置と速度の関係も同じ
x (t  t )  x(t )  v(t )t
t
 0
a( )d
積分
The Virtual Reality Society of Japan
位置・速度・加速度
x (t  t )  x(t )  v(t )t
x(0)  0
速度のときと同じように，位置を計算できる．
x(t)
v(t)
1
a(t)
2
3
The Virtual Reality Society of Japan
物理と微分積分
 物理学
 現実世界の性質（法則）を調べる学問
 現実世界の法則
 微分方程式で表される場合が多い
 法則の例
 運動方程式
f (t )  m
dv (t )
dt
2
2
d
h
(
x
,
t
)
d
h ( x, t )
2
波

c
dx2
dt 2
 コイルの電圧 V (t )  L di (t )
dt
h(x,t):波の高さ
c:波の伝わる早さ
微分方程式は簡単な式なのに，たくさんのことを表している．
The Virtual Reality Society of Japan
力学
 物理シミュレーション
 現実の物体は物理法則に従う．
 これをシミュレーションすれば，リアリティの高いVR世界が作れる．
 運動の法則
 運動方程式
力
f (t )  ma(t )
加速度
質量（比例係数）
 慣性の法則
dv (t )
⇔
a (t ) 
dt
f
a
v(t  t )  v(t )  a(t )t
 位置と速度の関係
x
(
t


t
)
dx (t )
⇔
v(t ) 
dt
加速度aが速度vを変える
 x(t )  v(t )t
速度vが位置xを変える
The Virtual Reality Society of Japan
運動のシミュレーション
f (t )  ma(t )
v(t  t )  v(t )  a(t )t
x(t  t )  x(t )  v(t )t
v(t  t )  v(t )  (f (t ) / m)t
x(t  t )  x(t )  v(t )t
f(t)が分かれば，車のときと同じようにシミュレーションできそう．
バネにつながった質点（おもり）の動きをシミュレーションしてみる
f(t)
f (t )  kx(t ) バネの法則
O
x
v(t  t )  v(t )  (kx(t ) / m)t
x(t  t )  x(t )  v(t )t
The Virtual Reality Society of Japan
運動のシミュレーション（1次元）
v(t  t )  v(t )  (kx(t ) / m)t
x(t  t )  x(t )  v(t )t
v(0)  0, x(0)  1 (最初の状態) k=1,m=1 t=0.1
v(0.1)  v(0)  x(0)0.1  0.1
0
1
v(0.2)  v(0.1)  x(0.1)0.1  ???
-0.1
???
vだけで計算することはできない．
x(0.1)  x(0)  v(0)0.1  1
v(0.2)  v(0.1)  x(0.1)0.1  0.2
-0.1
1
x(0.2)  x(0.1)  v(0.1)0.1  0.99
1
-0.1 xも求めれば計算できる．
The Virtual Reality Society of Japan
運動のシミュレーション（1次元）
シミュレーション結果
3
2
1
0
-1
0
5
10
15
20
速度v
位置ｘ
-2
-3
バネにつけたおもりの揺れがだんだん大きくなる．
そんなの見たことない．なにか変では？
モデルが変？誤差？ t=0.1のせい？
The Virtual Reality Society of Japan
運動のシミュレーション（1次元）
シミュレーション結果2 t=0.01
3
2
1
0
-1
0
5
10
15
20
速度v
位置ｘ
-2
-3
さっきより良くなった．誤差のせいだったようだ．
tは小さい方が良い．
でも無限に小さくすると無限に計算時間がかかる．
The Virtual Reality Society of Japan
運動のシミュレーション２（1次元）
fs(t)
O
f s (t )  kx(t )
f d (t )  bv(t )
x
バネの法則
ダンパーの法則
fd(t)
v(t  t )  v(t )  (kx(t )  bv(t ) / m)t
x(t  t )  x(t )  v(t )t
最初の状態： v(0)  0, x(0)  1
k  1, b  1
The Virtual Reality Society of Japan
運動のシミュレーション２（1次元）
シミュレーション結果 t=0.01
3
2
1
0
0
5
10
15
20
速度v
位置ｘ
-1
-2
-3
ダンパを入れると振幅がだんだん小さくなる．
The Virtual Reality Society of Japan
物理シミュレーション３（２次元）
形を持ったものをシミュレーション
おもりが壁にぶつかると
ぶつかっている間だけ，
バネ・ダンパを付ける．
p2 v2
m=1,k=1,b=1
1(自然長)
2
2
p3 v3
p1 v1
3
The Virtual Reality Society of Japan
物理シミュレーション３（２次元）
p2
自然長=1
バネによる力 f12s= -k s
(p1  p 2 )
s= (p1  p 2 ) 
| p1  p 2 |
e
y
O
ex
p1
v’2
p2
s のび=1

(p 1  p 2 )  (p 1  p 2 )


v

v’1 =  1
| p1  p 2 |  | p1  p 2 |

y
ex
px2  p y2
ダンパによる力 f12d= -b(v’1-v’2)
v2
e
O
| p | pの長さ 
v1
p1
v’1 f12
a  b  aと bの内積
| a || b | cos  a xbx  a y by

|b|cos
b
a
The Virtual Reality Society of Japan
物理シミュレーション３（２次元）
ey
O
f12  k ((p1  p 2 ) 
p2 v2
ex
p3 v3
p1v1
f13
f12
(p1  p 2 )
)
| p1  p 2 |

(p  p 2 )  (p1  p 2 )

 b ( v1  v 2 )  1
| p1  p 2 |  | p1  p 2 |

(p  p 3 )
f13   k ((p1  p 3 )  1
)
| p1  p 3 |

(p1  p 3 )  (p1  p 3 )


 b ( v 1  v 3 ) 
| p1  p 3 |  | p1  p 3 |

v1 (t  t )  v1 (t ) 
f12 (t )  f13 (t )
t
m
p1 (t  t )  p1 (t )  v1 (t )t
The Virtual Reality Society of Japan
物理シミュレーション３（２次元）
壁にぶつかっている場合：
2
p1 v1
d1R
 v1x 
f1R  kd1R  b 
0
f12 (t )  f13 (t )  f1R (t )
v1 (t  t )  v1 (t ) 
t
m
ほかの壁も同様
The Virtual Reality Society of Japan
物理シミュレーション３（２次元）
シミュレーション結果
ソースコード http://springhead.info/~hase/soft/sim2d.tar.gz
The Virtual Reality Society of Japan
剛体の運動
 剛体
 さっきの例のように，質点同士が動かないように固定されているもの
 バネダンパでもシミュレーションできるが，
剛体は剛体としてシミュレーションした方が効率が良い
 剛体の特徴
=
どんな移動も
並進（平行移動）と回転で表せる
並進と回転だけで運動を表そう
The Virtual Reality Society of Japan
剛体にかかる力
並進の力
運動方程式から
fii  fei  miai
f  f   m a
f  f   m a
f   m a
m1
p1 v1 a1
fi2
ii
fi1
p2 v2 a2 fe1 m5
m2
m3
m4
ii
fi2
i i
ei
ei
i i
i i
慣性の法則から
v(t  t )  v(t )  a(t )t
fe:外から加わった力（外力）
fi:質点同士が押し合う力（内力）
fi1
ei
f i1  fi 2
押したぶんだけ，押し返される
作用・反作用の法則
 m v (t  t )   m v (t )   m a (t )t
 m v (t  t )   m v (t )   f t
i
i
i
i
i
i
i
i
i i
ei
質量・速度 (運動量)の変化は，外から加えた力の和
The Virtual Reality Society of Japan
回転の力（トルク）（２次元）
r2
r1
O
f2
f1
回転軸（移動しない）
pr
fr
ex
f1のトルク：
f2のトルク：
r1 f1
r2 f2
fr = fのうち，質点をOを中心に回そうとする力
p
O
ey
力f1 f2は棒をまわそうとする．
まわす強さ：トルク
f
（残りは回転軸が力を出して打ち消す）
fのモーメント:|p||fr | = |pr||f | =p×f
p×f=
px fx
× =pxfy-pyfx
py fy
外積については布川昊, “２次元ベクトルの外積の効用(線形代数学の教科内容の改善に向けて) ”
http://www.dt.takuma-ct.ac.jp/~sawada/math/danwa5html/node14.html
を参照
The Virtual Reality Society of Japan
剛体にかかる回転の力（２次元）
内力のモーメントは打ち消しあう
fi1
運動方程式から
f ii  fei  mi ai
p i  (f ii  fei )  p i  mi ai
p1 p1r
p f  p f  p  m a
p f  p  m a
i
ii
i
ei
i
ei
i
i
p2
fi2
i i
0
i i
慣性の法則から
v i (t  t )  v i (t )  ai (t )t
O
p2r
p1  f i1 | p1r || f i1 |
  | p 2r || f i2 | p 2  fi 2
p i  v i (t  t )  p i  ( v i (t )  ai (t )t )
 p  m v (t  t )   p  m v (t )   p  m a (t )t
 p  m v (t  t )   p  m v (t )   p  f t
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i i
ei
位置×(質量・速度) (角運動量)の変化は，外から加えた位置×力(トルク)の和
The Virtual Reality Society of Japan
剛体の速度
速度
vg
＝
O
速度も並進と回転に分けられる vi =wp + vg
回転の中心Oはどこでも良い
vi
wp
The Virtual Reality Society of Japan
剛体の速度（２次元）
角速度：角度の速度
v

 (t  t )   (t )  (t )t
角速度

O
vg
O
pi
wi
vi
p P
v= 90°回転( p )
vi =wi + vg
vi= 90°回転( pi )+ vg
剛体の点piの速度は，
回転速度 と並進速度 vで表せる
The Virtual Reality Society of Japan
剛体の運動の式（並進）
m v (t  t)  m v (t)  f
i
i
i
i
ei
t
質量・速度 (運動量)の変化は，外から加えた力の和
剛体上の点piの速度vi = 90°回転( pi ) + vg
 mi v i   mi (90回転(pi )  v g )
  mi v g   90回転( mi p i )
( mipi  0となるように原点を取
＝ mi v g
M Mは剛体の質量と呼ばれます
Mv g (t  t )  Mv g (t )  fei t
剛体の運動の式（並進）
ると )
The Virtual Reality Society of Japan
剛体の運動の式（回転）（２次元）
p  m v (t  t)  p  m v (t)  p  f
i
i
i
i
i
i
i
ei
t
位置×(質量・速度) (角運動量)の変化は，外から加えた位置×力(トルク)の和
剛体上の点piの速度vi = 90°回転( pi ) + vg
 pi  mi v i   pi  mi (90回転(pi )  v g )
 ( mi p i )v g   mi p i  90回転(p i ) 
( mipi  0となるように原点を取
ると )
 ( mi | p i |2 )
I
Iは慣性モーメントと呼ばれます
Iω(t  t )  Iω(t )  pi  fei t
剛体の運動の式（回転）
The Virtual Reality Society of Japan
剛体の運動の式（２次元）
並進と回転をまとめると
Mv g (t  t )  Mv g (t )  fei t
Iω(t  t )  Iω(t )  pi  fei t
速度と位置の関係
x g (t  t )  x g (t )  v g (t )t
θ g (t  t )  θ g (t )  ωg (t )t
剛体の運動方程式
Mv g (t  t )  Mv g (t )   fei t
M
M
これで，シミュレーションできる
v g (t  t )  v g (t )
t
dtv g (t )
dt
  fei
  fei
Iω(t  t )  Iω(t )   p i  f ei t
Iω(t  t )  Iω(t )
  p i  f ei
t
dω(t )
I
  p i  f ei
dt
The Virtual Reality Society of Japan
3次元の剛体
2次元では，回転軸は面に垂直
だから，角速度・トルクは1次元

3次元では？
回転軸は3次元ベクトル
角速度も3次元ベクトル
 の向き：回転軸の向き． ||：回転速度
 v= ×p
トルクも3次元ベクトル
N =p×f
N
p
f
p
v
The Virtual Reality Society of Japan
3次元の外積
a×b
 3次元版 a×bとは？
|a×b|=|a||b’| = |a’||b| 2次元のときと同じ
a×bの向きは
aがx軸，b’がy軸のときのz軸の向き
b
a
数値で書くとこうなります．
 a x   bx   a y bz  by a z   0
 
    
a  b   a y    by    a z bx  bz a x    a z
a  b  a b b a   a
 z  z  x y x y  y
a  b  b  a
外積の計算については
 az
0
ax
a’ b’
a y  bx 
 
 a x  by 
0  bz 
布川昊, “２次元ベクトルの外積の効用(線形代数学の教科内容の改善に向けて) ”
http://www.dt.takuma-ct.ac.jp/~sawada/math/danwa5html/node14.html を参照
線形代数の本が外積について触れないのは，4次元以降が大変だから
ほんとに知りたければ http://members.jcom.home.ne.jp/1228180001/ をどうぞ．
The Virtual Reality Society of Japan
3次元の慣性モーメント

剛体上の点piの速度vi = ×pi + vg
 p  m v   p  m (ω  p  v
 ( m p )v   m p  (ω  p )
  m p  (ω  p )
i
i
i
i
i
i
i
i
g
i
i
i
i
i
g
pi
)
i
pi××pi 
 0
 p z p y  0
 pz p y 



  mi  p z
0  px  p z
0  p x ω
 p
  p

p
0
p
0
x
x
 y
 y

 p y 2  pz 2  px p y
 px pz 


2
2
  mi   p y p x
pz  px
 p y p z ω  Iω

2
2


p
p

p
p
p

p
z x
z y
x
y 

×pi

pi
×pi
I:慣性モーメント
The Virtual Reality Society of Japan
3次元の剛体の運動
m v (t  t)  m v (t)  f
i
i
i
i
ei
t
質量・速度 (運動量)の変化は，外から加えた力の和
Mv g (t  t )  Mv g (t )  fei t
p  m v (t  t)  p  m v (t)  p  f
i
i
i
i
i
i
i
ei
t
位置×(質量・速度) (角運動量)の変化は，外から加えた位置×力(トルク)の和
I(t  t )ω(t  t )  I(t )ω(t )  pi  fei t
速度と位置の関係
x g (t  t )  x g (t )  v g (t )t
θ g (t  t )  θ g (t )  ωg (t )t
The Virtual Reality Society of Japan
3次元の向きと角度
3次元での向きの表し方
回転ベクトル
y
y
x
z p/2,0,0
||:回転角， の向き：回転軸
足せない．
y
 p/2,0,0のあと 0,p/2,0
 p/2, p/2,0
基底(ex ey ez)
y
元の座標系で，その向きの
基底を数値にしたもの
変換の変換と考えれば足せる．
x
z
0,p/2,0
x
z
x
z
p/2, p/2,0
y
x
z
The Virtual Reality Society of Japan
回転のシミュレーション
 剛体の向きを表す基底(ex ey ez)を並べた行列Rを考える．

e
y
wy
ex
R
wx
ez
wz
t だけ回転する． 回転する行列D(dx dy dz)を用意する
ey
wy
Rの数値の基底をDだと考えると
D
ex RはDだけ回ることになる．

R
wx
R(t+t)= DR(t)
ez
w
z
The Virtual Reality Society of Japan
慣性モーメントの変換
I(t  t )ω(t  t )  I(t )ω(t )  pi  fei t 回転の運動
 Iは
 p y 2  pz 2

 mi   p y px
  pz px

 px p y
pz  px
 pz p y
2
2
 px pz 

 p y pz 
2
2
p x  p y 
なので，物体が回ると変化する．
 そこで，物体の基底Rでの慣性モーメントIRを考える．
I Rω R (t  t )  I Rω R (t )   pi R  fei R t
I R R 1ω(t  t )  I R R 1ω(t )  R 1 ( pi  fei t )
RI R R 1ω(t  t )  RI R R 1ω(t )   pi  fei t
I
I  RIR R 1
The Virtual Reality Society of Japan
3次元剛体運動のシミュレーション
並進
Mv g (t  t )  Mv g (t )  fei t
x g (t  t )  x g (t )  v g (t )t
v g (t  t ) を求める
x g (t  t ) を求める
回転
I(t  t )ω(t  t )  I(t )ω(t )  pi  fei t
I(t  t )ω(t  t ) を求める
I(t  t )  R(t )I R R(t )1 から，ω(t  t ) を求める
t だけ回転する． 回転する行列D(dx dy dz)を用意する
R(t+t)= DR(t)
The Virtual Reality Society of Japan
3次元剛体シミュレーションの例
ソースコード
http://springhead.info/
The Virtual Reality Society of Japan
バネダンパと制御
fs(t)
O
f s (t )  kx(t )
f d (t )  bv(t )
x
バネの法則
ダンパーの法則
fd(t)
v(t  t )  v(t )  (kx(t )  bv(t ) / m)t
x(t  t )  x(t )  v(t )t
最初の状態： v(0)  0, x(0)  1
k, b いろいろ
The Virtual Reality Society of Japan
いろいろなk,bでのシミュレーション
シミュレーション結果 t=0.1
2
2
1.5
1.5
1
1
0.5
0.5
0
-0.5 0
5
10
15
20
速度v
位置ｘ
0
-0.5 0
-1
-1
-1.5
-1.5
-2
-2
安定する
k=1,b=1
2
2
1.5
1
1
0.5
0
5
10
15
20
速度v
位置ｘ
15
20
0.5
0
-0.5 0
-1
-1
-1.5
-1.5
-2
-2
k=2,b=0.1
10
k=1,b=0.1
1.5
-0.5 0
5
速度v
位置ｘ
発散してしまう
5
10
15
k=4,b=0.2
20
速度v
位置ｘ
The Virtual Reality Society of Japan
安定する条件
 式を行列で書く
v(t  t )  v(t )  (kx(t )  bv(t ) / m)t
x(t  t )  x(t )  v(t )t
 v(t  t )  1  tb / m  tk / m  v(t ) 

  


t
1
 x(t  t )  
 x(t ) 
x(t  t )
x(t )
A
x(nt )  Ax((n  1)t )  AAx((n  2)t )
 A n x(0)
1  tb / m  tk / m 

 
t
1


n
 0
 
1
The Virtual Reality Society of Japan
安定する条件
 v(nt )  1  tb / m  tk / m 

  

t
1
 x(nt )  

n
0
 
1
An
Aという変換を次々にかけていくと，シミュレーションが進む．
n→∞のとき，Anがまともな値になればよい
b1
wy
a1
wx
A=(a1,b1)
Ai =(ai,bi)
と書くことにする．
The Virtual Reality Society of Japan
Anを考える
wy
A=(a1,b1)
Ai =(ai,bi)
b1
b2
b2
b1
a1
b3
a2
wx
a3
a1
a2
wy
重ねて書くと
b1
多分A∞は
a1
wx
･･･
00
00
安定だ
The Virtual Reality Society of Japan
Anを考える
Aの基底ベクトルが短ければAnは0になるのかな？
wy
b1
a1
a1もb1もwx wyよりは小さいが，
だんだん大きくなる．
wx
The Virtual Reality Society of Japan
Anを考える
wy
 1 0.5 

A  
 0.5 1 
b1
a1
wx
The Virtual Reality Society of Japan
Anを考える
wy
ためしに，
 1 0.5 

A  
 0.5 1 
b1
でrからArへの矢印を書いてみた
a1
wx
The Virtual Reality Society of Japan
Anを考える
wy
ためしに，
 1 0.5 

A  
 0.5 1 
b1
でrからArへの矢印を書いてみた
a1
wx
r
Ar矢印1個進む
A2rで矢印2個進む
∴ Anrは矢印の行き着く先
The Virtual Reality Society of Japan
Anを考える
wy
ためしに，
 1 0.5 

A  
 0.5 1 
b1
でrからArへの矢印を書いてみた
a1
wx
Ar矢印1個進む
A2rで矢印2個進む
∴ Anrは矢印の行き着く先
向きが変わらない場所がある
固有ベクトルと呼ぶ
The Virtual Reality Society of Japan
Anを考える
最初から，固有ベクトルが基底だと思ってみる
py
px
The Virtual Reality Society of Japan
Anを考える
最初から，固有ベクトルが基底だと思ってみる
実はただの拡大縮小
拡大率：1.5倍, 0.5倍
py
拡大率を固有値と呼ぶ
px
The Virtual Reality Society of Japan
Anを考える
 1 0.5 

A  
0
.
5
1


1
  1
固有ベクトルは p x   と p y   
1
1
1  1

P  
1
1


r2=Ar1 r1は(wx wy )を基底に数値にしてた． r’2 = 1.5 0 r’1
0 0.5
r’1 ：(px py )を基底にr1を数値にしたもの．
1.5 0 n
r2=Ar1
r’n = 0 0.5 r’1
r1 =Pr’1 r2=Pr’2
Pr’2 =APr’1
1.5n 0
r’2 =P-1APr’1
r’ =
r’
1
1  1  1 0.5 1  1
 


P AP  
1
1
0
.
5
1
1
1

 


 0.5 0.5  1 0.5 1  1



 

0
.
5
0
.
5
0
.
5
1
1
1




n
0 0.5n
1
1
 0.75 0.751  1 1.5 0 

  

 
  0.25 0.251 1   0 0.5 
1.5n 0
rn =Pr’n =P 0 0.5n r’1
n
=P 1.5 0n P-1 r1
0 0.5
固有値（拡大率）が
全部1未満なら安定
The Virtual Reality Society of Japan
固有値を見つけるには
 Ap=lpとなるようなlを見つければよい．
Ap=lp


1 0
 Eは単位行列 E  


のことです．
(A-lE)p=0



0 1


A-lE=0 ということは，
A-lEがぺちゃんこ変換⇔基底が行列の次元より少ない
ey
面積1
ex こんなのではなく，
面積1.5くらい
面積1
こんなの
面積0
ぺちゃんこ⇔面積0，これで書いたら計算できるかも
The Virtual Reality Society of Japan
行列式 determinant
 基底ベクトルが作る平行四辺形の面積をdetAと書いて，行列式とか
determinantと呼ぶ
 2次元の場合なら平行四辺形の面積・3次元なら平行6面体の体積
4次元なら平行8面体の体積．．．
ez
ey
ex
det(ex ey ) =ex×ey
ey
ex
det(ex ey ez)=(ex×ey )･ez
 ただし，detAは負にもなります
ey
ex
det(ex ey )は負
ex
ey
det(ex ey )は正
ey
ex
det(ex ey )は負
The Virtual Reality Society of Japan
行列式の計算
 4次元以上の場合は，
e2
 2 0.5 

det(e1 e 2 )  det
0
.
5
1


ke1
e1
e’2
wy
e’2
e2
wx
e1
 2 
 det  
  0.5 
e2
e3
のように変形し det(e1 e2) = e1x ･ e1yと計算できる
 2
2 6 2 4 
 


1
1 6 7 5 
e 4 )  det 

det
 3
3 10 6 10
 


2 7 6 9 
 2


 
 2 0  0
0

   
 
 1  3  6
 3
 det 

2
1
3 1  3

   
 
1
 2 1  4
   
 

 2  3  1  ( 2)  12
0 
 2
  1.75
 det
0
.
5
0
.
875


e1
 4次元の場合の計算例
det (e1
 0.5 
 2 
   0.25  
 1 
 0.5  
 0 0
2
   

 3  3
1

1

det

 4 1
3
   

5 1
2
   

6
 2
 
 
6
1

3
10
 3
 
 
7
 2
 
 
0
0
3
1
0
1
1
2
2
0


0
1

det
3
3


2
4 

 2  2
   
7 1
 6   1 3 
   
 6  2
   
0
3
1
1
0
 
0
1
 
 2
 
4
 2
2
 
 

5
1
1

2

det

10
 3
3
 
 

9
 2
2
 
 

0
0
2
 
 

0
0
1

3

det

 3
1
3
 
 

 4
 2
2
 
 

0
0
3
6
1
1
3
4
0
0
3
1
0
1
1
2
0

3
4

5 
0 

0 
0 

 2 
The Virtual Reality Society of Japan
固有値を見つけるには
Ap=lpとなるようなlを見つければよい．
Ap=lp (A-lE)p=0 A-lE=0 det(A-lE)=0
 1 0.5 

A  
 0.5 1  のとき，
 1 0.5   l 0 
  
)
det(A  lE)  det(
 0.5 1   0 l 
1  l 0.5  1  l   0.5 
  
  

 det
 0.5 1  l   0.5  1  l 
 (1  l ) 2  0.5 2  l2  2l  0.75  (l  1.5)(l  0.5)  0
 l  1.5, 0.5
The Virtual Reality Society of Japan
安定する条件
x(nt )  A n x(0)
1  tb / m  tk / m 

 
t
1


n
0
 
1
1  tb / m  tk / m 
1 0
1  tb / m  l
  l 
)  det
det(
t
1
t


0 1

 (1  tb / m  l )(1  l )  (tk / m)(t )
 l2  (2  tb / m)l  1  tb / m 
~ bt ~ kt 2
(b 
,k 
とおいた )
m
m
~
~ 2
~ ~
~
~
~
2  b  (2  b )  4(1  b  k ) 2  b  b 2  4k
l

2
2
~
~ ~
 l  ( 2  b )l  1  b  k
2
 tk / m 

1 l 
The Virtual Reality Society of Japan
安定する条件
~
~
~
2  b  b 2  4k
l
, | l | 1
2
~
~
~
~
~
b  2  b 2  4k
| l |
1
2
~
~
k
b  2
2
~
i) b 2  4k  0 のとき
~
~
~ ~
b 2  4k  i 4k  b 2
~
~ ~
( 2  b ) 2  4k  b 2
~ ~
| l |
 1 b  k  1
2
~ ~
b k
~
~
iii) b 2  4k  0 かつ 2  b  0のとき
~ ~
b k
~
b
~
~
k
b  2
2
4
~
ii) b 2  4k  0かつ 2  b  0のとき
~
~
~
2  b  b 2  4k
| l |
1
2
~
k 0
2
~
k
~ ~
4k  b 2
4
The Virtual Reality Society of Japan
試してみると
50
40
30
20
10
0
-10 0
-20
-30
-40
-50
50
40
30
20
10
0
-10 0
-20
-30
-40
-50
50
40
30
20
10
0
-10 0
-20
-30
-40
-50
0.1
0.2 0.3
0.5
0.4 0.5 0.6
1
0.7 0.8
1.5
0.9
1
2
速度v
位置ｘ
速度v
位置ｘ
速度v
0.5
1
1.5
2
位置ｘ
50
40
30
20
10
0
-10 0
-20
-30
-40
-50
50
40
30
20
10
0
-10 0
-20
-30
-40
-50
~
b
4
速度v
0.5
1
1.5
2
位置ｘ
2
速度v
0.5
1
1.5
2
位置ｘ
4
~ bt ~ kt 2
b
,k 
m
m
~
k
The Virtual Reality Society of Japan
まとめ
 位置の数値化
 座標系の変換と行列
 逆変換
 回転変換
 行列のランク
 物理法則と差分方程式
 運動方程式（質点・剛体）
 安定性と行列のn乗
ポヒマスセンサ
CG
ロボットアーム
シミュレーション
制御
をざっと見てきました．バス旅行で車窓から眺めたようなものです．
ぜひ，自分の足で歩いてみてください．
The Virtual Reality Society of Japan
明後日から
 自分の問題を解くには…
 実際の問題と教科書と紙と鉛筆を用意して考えてみてください．
→自分の足で歩く
 答えに自信がないときは，Excelか適当なプログラミング言語で確認
してみてください．
 お勧めの教科書
 線形代数：
平岡和幸, 堀玄(著) 「プログラミングのための線形代数」
 微分積分:
二見靖彦(著) 「理工系のための初等解析学とその応用」
 力学：
小出昭一郎(著) 「力学」物理テキストシリーズ
 3DCGと物理シミュレーションに必要なトピック：
Eric Lengyel (著), 狩野智英 (訳)
「ゲームプログラミングのための3Dグラフィックス数学」

Download Report