インタラクティブ技術に必要な 数学と物理 東京工業大学 長谷川晶一.

Slides:



Advertisements
Similar presentations
シミュレーション演習 G. 総合演習 ( Mathematica 演 習) システム創成情報工学科 テキスト作成: 藤尾 光彦 講義担当: 尾下 真樹.
Advertisements

Absolute Orientation. Absolute Orientation の問題 二つの座標系の間における剛体 (rigid body) 変換を復元す る問題である。 例えば: 2 台のステレオカメラから得られた3次元情報の間の関 係を推定する問題。 2 台のステレオカメラから得られた3次元情報の間の関.
1 線形代数学. 2 履修にあたって 電子情報システム学科 必修 2005 年度1セメスタ開講 担当 草苅良至 (電子情報システム学科) 教官室: G I 511 内線: 2095 質問等は上記のいずれかに行なうこと。 注意計算用のノートを準備すること。
1 運動方程式の例2:重力. 2 x 軸、 y 軸、 z 軸方向の単位ベクトル(長さ1)。 x y z O 基本ベクトルの復習 もし軸が動かない場合は、座標で書くと、 参考:動く電車の中で基本ベクトルを考える場合は、 基本ベクトルは時間の関数になるので、 時間で微分して0にならない場合がある。
陰関数定理と比較静学 モデルの連立方程式体系で表されるとき パラメータが変化したとき 如何に変数が変化するか 至るところに出てくる.
第1回 確率変数、確率分布 確率・統計Ⅰ ここです! 確率変数と確率分布 確率変数の同時分布、独立性 確率変数の平均 確率変数の分散
情報処理演習 (9)グラフィックス システム科学領域 日浦 慎作.
CGアニメーションの原理 基本技術 対象物体の動きや変形の設定方法 レンダリング技術
電磁気学C Electromagnetics C 7/27講義分 点電荷による電磁波の放射 山田 博仁.
・力のモーメント ・角運動量 ・力のモーメントと角運動量の関係
伝達事項 皆さんに数学と物理の全国統一テストを受けても らいましたが、この時の試験をまた受けていただ きます。
復習.
コリオリ力の復習資料 見延 庄士郎(海洋気候物理学研究室)
剛体の物理シミュレーション は難しい? 佐藤研助手 長谷川晶一.
スペクトル法による数値計算の原理 -一次元線形・非線形移流問題の場合-
平成23年8月 情報学群 岡田 守 このスライドは, 前川佳徳編著による「コンピュータグラフィックス」(オーム社)を基に作成されている.
大阪工業大学 情報科学部 情報システム学科 宇宙物理研究室 B 木村悠哉
次に 円筒座標系で、 速度ベクトルと加速度ベクトルを 求める.
工業力学 補足・復習スライド 第13回:偏心衝突,仕事 Industrial Mechanics.
Mathematica入門 数学を数式処理システムで 上智大学理工学部 大槻東巳 TA: 吉本行気,清水元気 2012年6月.
線形代数学 4.行列式 吉村 裕一.
多変数関数の積分(6/3~24) 重積分(2重積分) 第6章(§5は除く) 重積分の定義 「連続関数は積分可能」
伝達事項 試験は6/6 (土) 1限目の予定です。.
エンタテインメントのための 力覚提示 東京工業大学 精密工学研究所 長谷川晶一.
電気回路学Ⅱ エネルギーインテリジェンスコース 5セメ 山田 博仁.
ストークスの定理と、 渦度・循環の関係を 直感で理解する方法
(ラプラス変換の復習) 教科書には相当する章はない
10. 積分 積分・・確率モデルと動学モデルで使われる この章は計算方法の紹介 積分の定義から
3次元での回転表示について.
動力学(Dynamics) 運動方程式のまとめ 2008.6.17
シミュレーション演習 G. 総合演習 (Mathematica演習) システム創成情報工学科
ロボティクス・メカトロニクスの基礎 東京大学 生産技術研究所 倉爪 亮.
東京大学大学院情報理工学系研究科 Y.Sawa
電界中の電子の運動 シミュレータ作成 精密工学科プログラミング基礎 資料.
電磁気学C Electromagnetics C 5/28講義分 電磁波の反射と透過 山田 博仁.
OpenGLライブラリを用いた3次元フラクタルの描画
細胞の形と変形のための データ駆動型解析手法
Curriki原典
独立成分分析 5 アルゴリズムの安定性と効率 2007/10/24   名雪 勲.
カオス水車のシミュレーションと その現象解析
6. ラプラス変換.
数学教育・工学教育における 数式処理電卓の活用
メンバー 梶川知宏 加藤直人 ロッケンバッハ怜 指導教員 藤田俊明
電磁気学C Electromagnetics C 7/17講義分 点電荷による電磁波の放射 山田 博仁.
物理学Ⅰ - 第 11 回 - 前回のまとめ 回転軸の方向が変化しない運動 回転運動のエネルギーとその応用 剛体の回転運動の方程式
3次元での回転表示について.
物理学Ⅰ - 第 9 回 -.
知能システム論I(13) 行列の演算と応用(Matrix) 2008.7.8.
主成分分析 Principal Component Analysis PCA
動力学(Dynamics) 力と運動方程式 2008.6.10
変換されても変換されない頑固ベクトル どうしたら頑固になれるか 頑固なベクトルは何に使える?
ポリゴンメッシュ (2) - 変形と簡略化- 東京大学 精密工学専攻 大竹豊 資料および授業の情報は :
パターン認識特論 担当:和田 俊和 部屋 A513 主成分分析
電磁気学Ⅱ Electromagnetics Ⅱ 6/9講義分 電磁場の波動方程式 山田 博仁.
CGと形状モデリング 授業資料 1,2限: 大竹豊(東京大学) 3,4限: 俵 丈展(理化学研究所)
平面波 ・・・ 平面状に一様な電磁界が一群となって伝搬する波
速度ポテンシャルと 流線関数を ベクトルで理解する方法
資料 線型変換のイメージ 固有値、固有ベクトル 平賀譲(209研究室) 資料
4. システムの安定性.
逆運動学:手首自由度 運動学:速度、ャコビアン 2008.5.27
行列式 方程式の解 Cramerの公式 余因数展開.
宿題を提出し,宿題用解答用紙を 1人2枚まで必要に応じてとってください 配布物:ノート 2枚 (p.85~89), 小テスト用解答用紙 1枚
大阪工業大学 情報科学部 情報システム学科 学生番号 B02-014 伊藤 誠
【第六講義】非線形微分方程式.
目で見る一次変換 河合塾 数学科 生越茂樹 オゴセ シゲキ.
バネモデルの シミュレータ作成 精密工学科プログラミング基礎 資料.
電気回路学Ⅱ 通信工学コース 5セメ 山田 博仁.
力覚インタラクションのための 物理ベースモデリング
宿題を提出してください. 配布物:ノート 3枚 (p.49~60), 中間アンケート, 解答用紙 3枚 (1枚は小テスト,2枚は宿題用)
Cプログラミング演習 ニュートン法による方程式の求解.
Presentation transcript:

インタラクティブ技術に必要な 数学と物理 東京工業大学 長谷川晶一

自己紹介 長谷川晶一 VRをはじめて11年.学部2年からVRコンテストに参加 コンテストで忙しくて授業に出ない日々 教養の数学・物理は,ぎりぎりで単位取得 VRコンテスト→東工大 佐藤・小池研 (力覚インタフェースSPIDARの研究室) 最近は動力学シミュレーションも研究

VRシステムの中の数学と物理 VRシステム 提示装置 コンピュータ 数学 物理 計測装置 大画面 HMD 3DCGレンダリング ポヒマス 加速度センサ 数学 物理 数学 物理 コンピュータ 3DCGレンダリング シミュレーション 数学 物理

今日の目標 数学がVR役立つということを納得してもらいたい. VRシステムを作るために必要な最小限の数学と 物理の話をする. あとで,自分の問題を解くために,教科書のどこを見ればよいか分かるように. 今日の話は(あたりまえですが)不完全です.いつか,それに気づいてほしいです.

話の進行 位置の数値化 座標系の変換と行列 積分・微分と シミュレーション 物理シミュレーション 安定性と行列のn乗 トピック アプリケーション 教科書 位置の数値化 座標系の変換と行列 逆変換 回転変換 行列のランク 積分・微分と シミュレーション 物理シミュレーション 運動方程式(質点・剛体) 安定性と行列のn乗 位置・角度センサ 3DCG ロボット シミュレーション 制御 線形代数 微分積分 力学

位置と数値 ポヒマスセンサ という位置センサ トランスミッタからみた,レシーバ位置を計測 レシーバ トランスミッタ

位置と数値 数値で表さなければならない コンピュータ レシーバ トランスミッタ

位置と数値 レシーバの位置(P)は,トランスミッタ(O)から 矢印pだけ進んだところ 矢印pを数値で表したい. P p O

位置と数値 レシーバの位置(P)は,トランスミッタ(O)から 矢印pだけ進んだところ 矢印pを数値で表したい. 向きと長さの基準が必要:基底(基準の矢印) P ez ex ey p O 矢印p は,ex 3つ分+ ey 2つ分+ ez 1つ分 これを,p=3ex + 2ey + 1ez と書く

位置と数値 レシーバの位置(P)は,トランスミッタ(O)から 矢印p =3ex + 2ey + 1ezだけ進んだところ 原点O,基底(ex ey ez ) が決まると, 点Pは3つの数値 で表せる. (原点O,基底(ex ey ez ) )のことを座標系といいます. P ez ex ey p O 3 2 1

位置と数値 p=3ex + 2ey + 1ez p = ez ex ey 基底のおかげで 矢印が数値で表せるようになった 3 2 1 ベクトル 3 2 1 p = 基底は(ex ,ey ,ez ) 基底を忘れてはならないのだけど... 面倒なので省略します. ついでにいうと  基底ベクトル ex ey ez も  矢印だから数値で書けます. それにも,もちろん,基底(ex ,ey ,ez )が必要です. 1 ex = ey = ez =

座標系 手の位置の計測 トランスミッタを目の位置に置くと 目にくっついた座標系で,手の位置が数値になる. ez ex ey トランスミッタ O ph Ph ph = 2ex -1ey +3ez レシーバ 2 -1 3

座標系の変換 机から見た手の位置 ポヒマスが2個あると 2 1 pe = 2wx +1wy +1wz ey ex wz wx wy 2 -1 3 ph = 2ex -1ey +3ez pe ph ez 目にくっついた座標系eで 手の位置が数値になる q q = pe+ph :机から見た手の位置 机にくっついた座標系wでの 手の位置を数値にすると?

向きの計測 ポヒマスセンサは実は向きも測れる (ex ey ez )=( ) 位置pと同様に向き(ex ey ez )も数値で返す wy コンピュータ トランスミッタ wx レシーバ wz ez ex ey wy 0.5 0.9 1 -0.9 (ex ey ez )=( ) wx 基底は(wx wy wz ) ez ex wz

基底の変換(取替え) ex = 0.5 wx + 0 wy + 0.9 wz ey = 0 wx + 1 wy + 0 wz ポヒマスの返す数値から wz wx wy ey 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 ez だから,eをwで書き直せて, ex ph 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 2 1 -1.7 -1 3.3 0.3 4.3 机にくっついた座標系での 手の位置の数値 q=pe+ ph = + =

基底の変換と行列 簡単に書きたい x ax+byを(a b) と書くことにする. y ph = 2ex -1ey +3ez ex = 0.5 wx + 0.9 wz ey = wy ez =-0.9 wx + 0.5 wz なので, ph = (ex ey ez ) 0.5 0.9 1 -0.9 0.5 2 -1 3 = (wx wy wz ) (wx wy wz ) (wx wy wz ) 0.5 0 -0.9 0 1 0 0.9 0 0.5 2 -1 3 -1.7 -1 3.3 = (wx wy wz ) = (wx wy wz ) 基底を数値で表したもの→行列

基底の変換と行列のまとめ ey phは,2ex +1ey +3ez phは, wy 2 1 3 wx これをwで表したい ez ex wz 0.5 0.9 1 -0.9 0.5 (ex ey ez )は,( ) なので,phは, p 2 1 3 0.5 0 -0.9 0 1 0 0.9 0 0.5 2 1 3 -1.7 -1 3.3 ph= = 行列 行列とは基底を数値で書いて並べたもの

いろいろな変換 2DCGの変形 wy wx ey ex ey ey ex ex ex ey 線の描画→ 始点と終点を数値で指定 線の描画→ 始点と終点を数値で指定 基底:画面上のwx wy 基底を変えると変形する. wx wy 1 (ex ey )=( ) ex ey 1 (ex ey )=( ) ex ey ex ey 2 1 (ex ey )=( ) ex 1 (ex ey )=( ) ey 1 (ex ey )=( ) 0.5 (ex ey )=( ) 1 -1 (ex ey )=( )

回転変換 wy ey ex wx 回転行列 というのがあったような 1 回転変換→回転した基底での数値をもとの基底での数値に変換する.                             というのがあったような 回転変換→回転した基底での数値をもとの基底での数値に変換する. wx wy 1 q ex ey sinq cosq

逆変換 机の上の1台のトランスミッタで 目の座標系で表した手の位置を計測したい ポヒマスの出力: pe q ex ey ezの数値 基底は(wx wy wz ) 知りたいもの: phの数値 基底は(ex ey ez ) ez ex ey wz wx wy ph pe ph1 ph2 ph3 q O ph=q-pe なので, 基底(wx wy wz )でのph数値 は わかる ph1w ph2w ph3w

逆変換 (wx wy wz)を基底にしたex ey ezの数値は ポヒマスから得られるので, wz ex ph1w ph2w ph3w = wx もし,逆に(ex ey ez) を基底としたwx wy wzの数値がわかれば, ph1 ph2 ph3 wx wy wz ph1w ph2w ph3w ez = wz ex ph1 ph2 ph3 なので, の数値がわかる. wx

逆変換 ph1w ph2w ph3w ex ey ez ph1 ph2 ph3 ph1 ph2 ph3 wx wy wz ph1w ph2w = = 逆行列 ph1w ph2w ph3w ex ey ez wx wy wz ph1w ph2w ph3w = ph1w ph2w ph3w 1 0 0 0 1 0 0 0 1 ph1w ph2w ph3w = ex ey ez wx wy wz 1 0 0 0 1 0 0 0 1 = -1 = ex ey ez wx wy wz と書く

逆行列の求め方 ex ey ez wx wy wz x2倍 z=z- 3/2 x z=z- 3/2 x 1/2 0 - 3/2 0 1 0 1/2 0 - 3/2 0 1 0 3/2 0 1/2 √ 1 0 0 0 1 0 0 0 1 √ x2倍 √ x2倍 1 0 - 3 0 1 0 3/2 0 1/2 2 0 0 0 1 0 0 0 1 √ z=z- 3/2 x √ √ 1 0 - 3 0 1 0 0 0 2 2 0 0 0 1 0 - 3 0 1 z=z- 3/2 x √ √

逆行列の求め方 z=z- 3/2 x z=z- 3/2 x ex ey ez wx wy wz √ 1 0 - 3 0 1 0 0 0 2 1 0 - 3 0 1 0 0 0 2 2 0 0 0 1 0 - 3 0 1 z=z- 3/2 x √ z=z- 3/2 x √ √ √ 1 0 - 3 0 1 0 0 0 1 2 0 0 0 1 0 - 3/2 0 1/2 1/2倍 √ z 1/2倍 x=x+ 3z √ √ 1 0 0 0 1 0 0 0 1 1/2 0 3/2 0 1 0 - 3/2 0 1/2 x=x+ 3z √ √ ex ey ez wx wy wz

3DCGでの変換の利用 ポリゴン・メッシュとは,3角形の頂点座標 W: ( )とO ワールド座標系 M:(mx my mz )とmp モデリングするときの座標系 V:(vx vy vz )とvp 視点座標系(視点の位置と向き) 1 p2 mz mx my p1 p3 vz vx vy 視点 wz wx wy mp vp O

3DCGの変換 p2 p3 p1 p1w= mp + Mp1m p1v =V-1 (mp +Mp1m-vp ) p1w=vp + Vp1v モデリングしたときの頂点位置p1 p2 p3は,M座標系での数値になっている.数値:p1m p2m p3m これをV座標系での数値:p1v p2v p3v にできれば画面に表示できる. p2 mz mx my p1 p3 vz vx vy 視点 wz wx wy mp vp O vz vx vy p1 p2 p3 p1w= mp + Mp1m p1v =V-1 (mp +Mp1m-vp ) p1w=vp + Vp1v p1v = V-1 (p1w -vp) p1v x p1v y p2v x p2v y あとは, , …に線を引けばよい

逆変換のない場合 戻せる ey ex 戻せない つぶれてしまうと元に戻せない 0.5 (ex ey )=( ) 1 -1 ex ey ex (ex ey )=( ) ex ey 戻せる 1 (ex ey )=( ) ex ey ex ey 2 1 (ex ey )=( ) ex ex ey ey (ex ey )=( ) 1 戻せない 0.5 (ex ey )=( ) 1 -1 つぶれてしまうと元に戻せない

行列のランク ロボットの関節の速度と手先の速度 e1: w1が1のときの手先の速度 v e1 r2 e2: w2が1のときの手先の速度 e2 wx wy O ロボットなので, r1, r2, q1, q2,はわかる. w1 , w2が設定できる. 手先に速度vを出したい r1

行列のランク w1 , w2がわかっているとき,手先の速度は となる. w→vができたのだから,逆変換 v → wをすれば, vを実現するwがわかる 逆変換できない場合は?

行列のランク e2 e1 基底が2つある. 逆変換あり. どんなvでも出せる. e1 e2 基底が1つしかない. 逆変換なし.vは直線上だけ r1=1, r2=1, q1=30°, q2 =90°のとき 基底が2つある. 逆変換あり. どんなvでも出せる. e1 r1=1, r2=0, q1=90°, q2 =90°のとき e2 基底が1つしかない. 逆変換なし.vは直線上だけ e2 e1 r1=1, r2=1, q1=30°, q2 =30°のとき 実は基底が1つしかない. 逆変換なし.vは直線上だけ

行列のランク e1 e2 基底がない 0個. 逆変換なし.vはまったく出せない 基底の数のことを行列のランクといいます. r1=0, r2=0, q1=30°, q2 =90°のとき e1 e2 基底がない 0個. 逆変換なし.vはまったく出せない 基底の数のことを行列のランクといいます. ・平行な基底は数えません. ・長さが0の基底は数えません. 2次元なら2つ,3次元なら3つ,4次元なら4つ …. 基底がないと,つぶれてしまい,逆変換できません

冗長マニピュレータと長方形行列 e1 e2 wy e3 O wx 速度vは2次元なのに,wは3つ 基底は3つあるように見えるが, e1 0.5 0.8 e2 = 0.5e1 +0.8e3 のように,ほかの2つで表せる 関節2を動かさなくても,関節1と3を動かせば, どんな速度vも出せる. 逆行列はない.答えが無数にある.

シミュレーション シミュレーション 例:車の動きのシミュレーション アクセル:原因 速度:結果 シミュレーション 速度 アクセルの踏み具合 現実のものの動きを計算機の中に再現すること そのために,現実世界の法則を見つけ,それを再現する. 例:車の動きのシミュレーション 「アクセルを踏むと加速する」→法則 アクセル:原因 速度:結果 シミュレーション 速度 アクセルの踏み具合 ? 時刻t=0で, 0m/sだった 時間t

「アクセルを踏むと車は加速する」法則を 式で書くと シミュレーション アクセル:原因 車の動き:結果 シミュレーション 「アクセルを踏むと車は加速する」法則を 式で書くと ちょっと後(t+Dt)の速度= 今(t)の速度+ 今のアクセル×踏んだ時間 今のアクセル t t +Dt Dtが小さければ a(t) と a(t+Dt)はだいたい等しい Dtが小さければ小さいほど誤差が少ない.

数と関数 t a0 a1 a2 a3 a4 a(t) ここまでは,aといったら変数だったので値は1つだった. a0 a1 a2 a3 a4 a(t) ここまでは,aといったら変数だったので値は1つだった. a(t)は,a0 a1a2 …を含んでいる. だけど,これからは,a(t)と書かないでaと書く.こともある. この場合,aは関数.たくさんの値を持っていて,tごとに 値a(t)がある.

これをプログラムでやるのが,コンピュータシミュレーション 法則 0.1 0.05 0.15 0.2 ... 0.4 0.45 分かっていること 最初の状態 2 1 その後の 経過 これをプログラムでやるのが,コンピュータシミュレーション

シミュレーションと積分 誤差をなくすためDtを0に近づける 積分

シミュレーションと積分 シミュレーションではDtずつ時間を進めて ほしい時刻まで繰り返し計算する. 積分はその結果のこと

シミュレーションと微分 車の速度v(t)の記録がある. アクセルの踏み具合a(t)を知りたい v(t) 0.275 0.225 1.0 1.1 v(t)がわかっているので, ちょっと後の値と今の値を引き算すれば出てくる.

シミュレーションと微分 誤差をなくすため,Dtを0に近づける 微分 vをtで微分する. vのt微分 微分が入った方程式 →微分方程式と呼ぶ

微分と積分 (ただしDtは十分小さい) これが元の式. a(t)とv(t)の関係を表している. 微分 積分 位置と速度の関係も同じ

位置・速度・加速度 速度のときと同じように,位置を計算できる. x(t) 2 3 1 v(t) a(t)

物理と微分積分 物理学 現実世界の法則 法則の例 微分方程式は簡単な式なのに,たくさんのことを表している. 現実世界の性質(法則)を調べる学問 現実世界の法則 微分方程式で表される場合が多い 法則の例 運動方程式 波 コイルの電圧 h(x,t):波の高さ c:波の伝わる早さ 微分方程式は簡単な式なのに,たくさんのことを表している.

力学 物理シミュレーション 運動の法則 f a 加速度 力 質量(比例係数) ⇔ 加速度aが速度vを変える ⇔ 速度vが位置xを変える 現実の物体は 物理法則 に従う. これをシミュレーションすれば,リアリティの高いVR世界が作れる. 運動の法則 運動方程式 慣性の法則 位置と速度の関係 f a 加速度 力 質量(比例係数) ⇔ 加速度aが速度vを変える ⇔ 速度vが位置xを変える

f(t)が分かれば,車のときと同じようにシミュレーションできそう. 運動のシミュレーション f(t)が分かれば,車のときと同じようにシミュレーションできそう. バネにつながった質点(おもり)の動きをシミュレーションしてみる f(t) バネの法則 O x

運動のシミュレーション(1次元) k=1,m=1 Dt=0.1 1 -0.1 ??? vだけで計算することはできない. -0.1 1 1 (最初の状態) k=1,m=1 Dt=0.1 1 -0.1 ??? vだけで計算することはできない. -0.1 1 1 -0.1 xも求めれば計算できる.

運動のシミュレーション(1次元) シミュレーション結果 バネにつけたおもりの揺れがだんだん大きくなる. そんなの見たことない.なにか変では? モデルが変?誤差? Dt=0.1のせい?

運動のシミュレーション(1次元) シミュレーション結果2 Dt=0.01 さっきより良くなった.誤差のせいだったようだ. Dtは小さい方が良い. でも無限に小さくすると無限に計算時間がかかる.

運動のシミュレーション2 (1次元) fs(t) バネの法則 O ダンパーの法則 x fd(t) 最初の状態:

運動のシミュレーション2 (1次元) シミュレーション結果 Dt=0.01 ダンパを入れると振幅がだんだん小さくなる.

物理シミュレーション3 (2次元) 形を持ったものをシミュレーション おもりが壁にぶつかると ぶつかっている間だけ, バネ・ダンパを付ける. p2 v2 m=1,k=1,b=1 1(自然長) 2 2 p3 v3 p1 v1 3

物理シミュレーション3 (2次元) p2 バネによる力 f12s= -k s ey s= O ex s p1 自然長=1 のび=1 バネによる力 f12s= -k s s= ダンパによる力 f12d= -b(v’1-v’2) v2 v’2 p2 v’1 = O ex ey v1 a b q p1 v’1 f12 |b|cosq

物理シミュレーション3 (2次元) O ex ey f13 p2 v2 p3 v3 p1v1 f12

物理シミュレーション3 (2次元) 壁にぶつかっている場合: 2 p1 v1 d1R ほかの壁も同様

物理シミュレーション3 (2次元) シミュレーション結果 ソースコード http://springhead.info/~hase/soft/sim2d.tar.gz

剛体の運動 剛体 剛体の特徴 = どんな移動も 並進(平行移動)と回転 で表せる 並進と回転だけで運動を表そう さっきの例のように,質点同士が動かないように固定されているもの バネダンパでもシミュレーションできるが, 剛体は剛体としてシミュレーションした方が効率が良い 剛体の特徴 = どんな移動も 並進(平行移動)と回転 で表せる 並進と回転だけで運動を表そう

剛体にかかる力 並進の力 m1 p1 v1 a1 m2 p2 v2 a2 m3 m4 m5 fi1 fi2 fe1 運動方程式から 慣性の法則から fe:外から加わった力(外力) fi:質点同士が押し合う力(内力) 押したぶんだけ,押し返される 作用・反作用の法則 fi1 fi2 質量・速度 (運動量)の変化は,外から加えた力の和

回転の力(トルク) (2次元) r2 r1 O 力f1 f2は棒をまわそうとする. まわす強さ:トルク f2 f1 回転軸(移動しない) pr fr = fのうち,質点をOを中心に回そうとする力 (残りは回転軸が力を出して打ち消す) fのモーメント:|p||fr | = |pr||f | =p×f p O f fr ex ey p×f= × =pxfy-pyfx px py fx fy 外積については 布川昊, “2次元ベクトルの外積の効用(線形代数学の教科内容の改善に向けて) ” http://www.dt.takuma-ct.ac.jp/~sawada/math/danwa5html/node14.html を参照

剛体にかかる回転の力(2次元) O 内力のモーメントは打ち消しあう 運動方程式から fi1 p1r p1 fi2 p2 p2r 慣性の法則から 位置×(質量・速度) (角運動量)の変化は,外から加えた位置×力(トルク)の和

剛体の速度 速度 = O vi wp vg 速度も並進と回転に分けられる 回転の中心Oはどこでも良い vi =wp + vg

剛体の速度(2次元) 角速度:角度の速度 q w O P v p v= 90°回転(w p ) pi vi vg vi =wi + vg vi= 90°回転(w pi )+ vg 剛体の点piの速度は, 回転速度w と並進速度 vで表せる

剛体の運動の式(並進) 剛体上の点piの速度vi = 90°回転(w pi ) + vg 剛体の運動の式(並進) M 質量・速度 (運動量)の変化は,外から加えた力の和 剛体上の点piの速度vi = 90°回転(w pi ) + vg M Mは剛体の質量と呼ばれます 剛体の運動の式(並進)

剛体の運動の式(回転) (2次元) 剛体上の点piの速度vi = 90°回転(w pi ) + vg 剛体の運動の式(回転) I 位置×(質量・速度) (角運動量)の変化は,外から加えた位置×力(トルク)の和 剛体上の点piの速度vi = 90°回転(w pi ) + vg I Iは慣性モーメントと呼ばれます 剛体の運動の式(回転)

剛体の運動の式(2次元) 並進と回転をまとめると 剛体の運動方程式 速度と位置の関係 これで,シミュレーションできる

3次元の剛体 2次元では,回転軸は面に垂直 3次元では? だから,角速度・トルクは1次元 回転軸は3次元ベクトル 角速度も3次元ベクトル w wの向き:回転軸の向き. |w|:回転速度 v= w×p トルクも3次元ベクトル N =p×f w p v N p f

3次元の外積 a×b 3次元版 a×bとは? |a×b|=|a||b’| = |a’||b| 2次元のときと同じ a×bの向きは aがx軸,b’がy軸のときのz軸の向き 数値で書くとこうなります. 外積の計算については 布川昊, “2次元ベクトルの外積の効用(線形代数学の教科内容の改善に向けて) ” http://www.dt.takuma-ct.ac.jp/~sawada/math/danwa5html/node14.html を参照 線形代数の本が外積について触れないのは,4次元以降が大変だから ほんとに知りたければ http://members.jcom.home.ne.jp/1228180001/ をどうぞ.

3次元の慣性モーメント w 剛体上の点piの速度vi = w×pi + vg pi w×pi w pi×(w×pi ) pi w×pi

3次元の剛体の運動 速度と位置の関係 質量・速度 (運動量)の変化は,外から加えた力の和 位置×(質量・速度) (角運動量)の変化は,外から加えた位置×力(トルク)の和 速度と位置の関係

3次元の向きと角度 3次元での向きの表し方 回転ベクトルq 基底(ex ey ez) |q|:回転角, qの向き:回転軸 足せない. (p/2,0,0)のあと (0,p/2,0) (p/2, p/2,0) 基底(ex ey ez) 元の座標系で,その向きの 基底を数値にしたもの 変換の変換と考えれば足せる. x x z z (p/2,0,0) (p/2, p/2,0) y x x y z (0,p/2,0) x z z y

回転のシミュレーション 剛体の向きを表す基底(ex ey ez)を並べた行列Rを考える. w ey wy ex R wx ez wz Dq=wDt だけ回転する. Dq回転する行列D(dx dy dz)を用意する wz wx wy ez ex ey Rの数値の基底をDだと考えると RはDだけ回ることになる. R(t+Dt)= DR(t) D Dq R

慣性モーメントの変換 回転の運動 Iは なので,物体が回ると変化する. そこで,物体の基底Rでの慣性モーメントIRを考える. I

3次元剛体運動のシミュレーション 並進 回転 R(t+Dt)= DR(t) を求める を求める を求める から, を求める Dq=wDt だけ回転する. Dq回転する行列D(dx dy dz)を用意する R(t+Dt)= DR(t)

3次元剛体シミュレーションの例 ソースコード http://springhead.info/

バネダンパと制御 fs(t) バネの法則 O ダンパーの法則 x fd(t) 最初の状態:

いろいろなk,bでのシミュレーション シミュレーション結果 Dt=0.1 安定する 発散してしまう k=1,b=1 k=1,b=0.1

安定する条件 式を行列で書く

安定する条件 An Aという変換を次々にかけていくと,シミュレーションが進む. n→∞のとき,Anがまともな値になればよい wy b1 a1 A=(a1,b1) Ai =(ai,bi) と書くことにする. wx

Anを考える wy a2 b2 a1 b1 b2 b3 b1 ・・・ a1 a3 wx a2 wy a1 b1 wx 0 0 安定だ A=(a1,b1) Ai =(ai,bi) b3 b1 ・・・ a1 a3 wx a2 wy a1 b1 wx 重ねて書くと 多分A∞は 0 0 安定だ

Anを考える Aの基底ベクトルが短ければAnは0になる のかな? wy b1 a1もb1もwx wyよりは小さいが, だんだん大きくなる.

Anを考える wy b1 a1 wx

Anを考える ためしに, wy b1 a1 でrからArへの矢印を書いてみた wx

Anを考える wy b1 a1 wx ためしに, でrからArへの矢印を書いてみた Ar矢印1個進む A2rで矢印2個進む ∴ Anrは矢印の行き着く先 r

Anを考える wy b1 a1 wx ためしに, でrからArへの矢印を書いてみた Ar矢印1個進む A2rで矢印2個進む ∴ Anrは矢印の行き着く先 向きが変わらない場所がある 固有ベクトルと呼ぶ

Anを考える 最初から,固有ベクトルが基底だと思ってみる a1 b1 wx wy py px

Anを考える py wy b1 px a1 wx 最初から,固有ベクトルが基底だと思ってみる 実はただの拡大縮小 拡大率:1.5倍, 0.5倍 拡大率を固有値と呼ぶ px

Anを考える r2=Ar1 r’1 :(px py )を基底にr1を数値にしたもの. r2=Ar1 Pr’2 =APr’1 r2=Pr’2 固有ベクトルは r2=Ar1 r1は(wx wy )を基底に数値にしてた. r’2 = r’1 1.5 0 0 0.5 r’1 :(px py )を基底にr1を数値にしたもの. n r2=Ar1 Pr’2 =APr’1 r’2 =P-1APr’1 1.5 0 0 0.5 r’n = r’1 r1 =Pr’1 r2=Pr’2 1.5n 0 0 0.5n r’n = r’1 1.5n 0 0 0.5n rn =Pr’n =P r’1 =P P-1 r1 1.5n 0 0 0.5n 固有値(拡大率)が 全部1未満なら安定

固有値を見つけるには Ap=lpとなるようなlを見つければよい. Ap=lp (A-lE)p=0 A-lE=0 ということは, ex ey こんなのではなく, こんなの 面積1 面積1.5くらい 面積1 面積0 ぺちゃんこ⇔面積0,これで書いたら計算できるかも

行列式 determinant ex ey det(ex ey ) =ex×ey det(ex ey ez)=(ex×ey )・ez ez 基底ベクトルが作る平行四辺形の面積をdetAと書いて,行列式とかdeterminantと呼ぶ 2次元の場合なら平行四辺形の面積・3次元なら平行6面体の体積 4次元なら平行8面体の体積... ただし,detAは負にもなります ex ey det(ex ey ) =ex×ey det(ex ey ez)=(ex×ey )・ez ez ey ex ex ex ey ey det(ex ey )は負 det(ex ey )は正 det(ex ey )は負

行列式の計算 4次元以上の場合は, e2 e1 4次元の場合の計算例 e’2 ke1 のように変形し det(e1 e2) = e1x ・ e1yと計算できる e2 e1 e’2 ke1 wx wy

固有値を見つけるには Ap=lpとなるようなlを見つければよい. Ap=lp ⇔(A-lE)p=0 ⇔ A-lE=0 ⇔ det(A-lE)=0 のとき,

安定する条件

安定する条件 iii) i) 4 ii) 2 4

試してみると 2 4

まとめ 位置の数値化 座標系の変換と行列 物理法則と差分方程式 安定性と行列のn乗 ポヒマスセンサ 逆変換 回転変換 CG 行列のランク 運動方程式(質点・剛体) 安定性と行列のn乗 ポヒマスセンサ CG ロボットアーム シミュレーション 制御 をざっと見てきました.バス旅行で車窓から眺めたようなものです. ぜひ,自分の足で歩いてみてください.

明後日から 自分の問題を解くには… お勧めの教科書 実際の問題と教科書と紙と鉛筆を用意して考えてみてください. →自分の足で歩く 答えに自信がないときは,Excelか適当なプログラミング言語で確認してみてください. お勧めの教科書 線形代数: 平岡 和幸, 堀玄(著) 「プログラミングのための線形代数」 微分積分: 二見靖彦(著) 「理工系のための初等解析学とその応用」 力学: 小出昭一郎(著) 「力学」物理テキストシリーズ 3DCGと物理シミュレーションに必要なトピック: Eric Lengyel (著), 狩野智英 (訳) 「ゲームプログラミングのための3Dグラフィックス数学」