[Unity] オイラー陽解法から始める物理シミュ入門
はじめに
Unity始めたばかり!という学生さんにも分かるように丁寧に書くよう心掛けました。
前提として 6. プログラムを実装 では実際に物理シミュレーションを試すためにUnityを使ったプログラミングの基本的な知識が必要です。しかし、前半はUnityが分からなくても読めるのではないかと思います。 なお、バージョンは Unity2021.3.14f1 を使っています。記事の内容ではバージョン依存はあまりないと思います。
(間違いなどあればご指摘いただけると助かります)
1. 概要
- 一番簡単で基本らしい 「オイラー陽解法」 をUnityで実装してみる
- 意味も用法も分からなかった教科書のテイラー展開の実用方法が分かる
- 「得たい物理の動き」の運動方程式が分かれば、ランタイムでシミュレーション計算して動かせる
いきなり運動方程式と言われても頭がフリーズする人は私を含めて結構多いと思います。
記憶が定かではないですが、高校の物理でやった
「得たい物理の動き」の運動方程式が分かれば
と書いた通り、運動方程式は既に世にある式を用いるので、あまり気負わないでください。あくまで 「道具として物理シミュレーションをする方法」 を説明したいと思います。
しかし、ただ道具としての使い方を説明するだけでなく、その際に用いられるテイラー展開などがどう使われるのか、その数式としての意味について解説し、物理シミュレーションを分かりやすく書ければと思います。
例として実装する振り子
2. 題材とする物理
運動方程式を用いて、その運動の物理シミュレーションをしていきます。
ではまず題材となりそうな運動方程式を探してきます。
振り子の運動方程式 は以下です。これを用いて物理シミュレーションをしていきましょう。
ただ、すぐには使わないので振り子の物理シミュレーションをするとだけ念頭において、ひとまず式は忘れて大丈夫です。
3. 少し未来の値を計算する
天気予報だったり、弾道予測だったり、シミュレーションによって未来の値を予測しますね。
今回は、現在時間
なお、
さて、未来予知の力を会得してないけど、どうやって未来
昔の数学の天才たちが、未来は分からないけれど特定条件下であれば未来の値が 近似 によって得られる方法を編み出してくれています。
冒頭で名前を出していますが、それが天才テイラーさんの名を冠した テイラー展開 です。
4. テイラー展開
関数はそのテイラー級数の有限個の項を用いて近似することができる。
点 a を含む実数の開区間 I ⊆ R 上で無限階微分可能な関数 f ∈ C∞(I) が与えられたとき、べき級数
\sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{k!}(x-a)^{k}
を関数 f の点 a まわりのテイラー級数という。
...どういうこと?
f(x) って?
4.1 そもそも関数 表記について、少しおさらいをしておきましょう。
仮に
プログラムの関数と同じですね。
int f(int x)
{
return x + 1;
}
4.2 テイラー級数って?
もし数式として関数
ですが、複雑な現実の物理においてリアルタイムで速度や力が変わってくる場合には単一の数式で計算することが難しい場合があります。
つまりは 「万能な
ですが、
関数はそのテイラー級数の有限個の項を用いて近似することができる。
そして、そのそれっぽい
魔法のように聞こえますが、まったく未知の
ちなみに上記の式を「テイラー級数」と呼びますが、「級数」というのは単に
例は以下です。「
このように、 欲しい
4.3 すごい式は分かったが...分からない
もう頭がパンクしそうですが、一言で言えば 「
とはいえ、結局はその数式がよく分からない...。ゆっくり見ていきましょう。
なお、テイラー級数は無限までたくさん足せばとても精度が高くなります。しかし、リアルタイムで計算するには現実的ではないので、今回は精度を犠牲に簡単さを取り、第2項までのとても簡素なものを扱います。
急にスリムになって寂しくなりました。でも見やすくなりましたね。続きを見ていきましょう。
4.4 ちょっと分かりやすい式になった...分かるかも
この
ほぼほぼ
ただし、テイラー級数の項が多ければその限りではありません。Wikipedia(テイラー展開)で項の数によって、どのぐらいの範囲で、どのぐらい精度が変わるか見れるのでぜひ見てみてください。
以上がテイラー展開の
4.5 やっと未来の値を手に入れる
式中の
では、これで式中の
とても分かりやすい式になりましたね!!
つまり、
簡単な数式
4.6 テイラー展開のまとめ
- 少し先の未来の値を近似式によって計算できる
-
自体は分からなくても、ある時点f(x) における値a とその微分f(a) が分かればいいf'(a)
5. 振り子の動きをシミュレーションしよう
やっと本題です。長かったですが、少し先の未来の値をテイラー展開を使って計算する雰囲気が掴めたでしょうか。近似式の
5.1 振り子の座標を計算する
続いて、 振り子の動きを思い出しておきましょう。
振り子の座標を動かしたいので、座標について考えます。
ですね。座標の時間微分
そもそも微分とはなんだったでしょうか。(Wikipedia(微分))
下図のように、
つまり、座標の時間微分とは 「極短い時間
よって、式は次のようになります。なお、時刻tにおける速度はVelocityのVをとって
...これは。見慣れている方もいそうです。速度が一定のときの等速直線運動をプログラムで書くときと同じですね。テイラー展開から始まったため難しそうに感じましたが普段から書いている座標計算です。
ただ、振り子の速度が分からない上に、振り子は等速直線運動ではないのでまだ計算することはできません。
しかし、あとは速度さえ分かれば計算ができます。
5.2 振り子の速度を計算する
ここでも同様にテイラー展開に頼ります。
現在の速度
速度の時間微分
そう、加速度です!!また単位を見てみましょう。加速度の単位は
となり、速度の時間微分も
よって、加速度さえわかればテイラー展開を使って未来の速度が計算できます。
5.3 振り子の加速度を計算する
ここでもまたテイラー展開...だと終わりが見えないですね。もう使わずとも実は手に入ります。 振り子の運動方程式 を思い出してみましょう。
ここでの
- 現在の座標:
(座標と言っていますが運動方程式にある円弧の長さのこと)x - 重力加速度の大きさ:
g - 糸の長さ:
l
これだけで加速度が求められます。そう、最初に出てきた運動方程式
これでやっと振り子の座標がシミュレーションで計算できそうです。
6. プログラムを実装
6.1 ひな形
順番にコードに落とし込んでいきましょう。まずはひな形を用意しておきましょう。
計算は Time.deltaTime
が固定で十分小さいことが保証されている FixedUpdate()
で行います。 Update()
で Time.deltaTime
を取ると、テイラー展開の条件だった Time.deltaTime
の値は
Update()
で計算することもできますが、経過時間分だけ固定 while
を使って1frameで数回のシミュレーションを回す必要があります。コードが増えるので、今回は素直に FixedUpdate()
を使います。
public class Pendulum : MonoBehaviour
{
void Start()
{
}
void FixedUpdate()
{
}
}
6.2 加速度を計算する
加速度、速度、座標の順で計算していけばよいですね。では、まずは加速度
加速度は
でした。
public class Pendulum : MonoBehaviour
{
// t+dt秒の計算をするためにt秒の情報が必要なので、取っておくための変数
float _x;
void Start()
{
// 初期位置をいれておく。開始時の振り子の円弧の長さ。
_x = 1;
}
void FixedUpdate()
{
// [説明変数] 定数
float dt = Time.deltaTime;
float g = Physics.gravity.magnitude;
float l = 5;
// [説明変数] 現在の値
float x_t = _x;
// 加速度(速度の時間微分)
float a = -g / l * x_t;
}
}
x_t
は更新されておらず初期値のままなので、加速度も固定のままです。しかし、最終的に x_t
が更新され始めれば、 加速度
加速度
6.3 速度を計算する
加速度
public class Pendulum : MonoBehaviour
{
// t+dt秒の計算をするためにt秒の情報が必要なので、取っておくための変数
float _x;
float _v; // <===== NEW
void Start()
{
// 初期位置をいれておく。開始時の振り子の円弧の長さ。
_x = 1;
}
void FixedUpdate()
{
// [説明変数] 定数
float dt = Time.deltaTime;
float g = Physics.gravity.magnitude;
float l = 5;
// [説明変数] 現在の値
float x_t = _x;
float v_t = _v; // <===== NEW
// 加速度(速度の時間微分)
float a = -g / l * x_t;
// dt秒先の速度を計算
_v = v_t + a * dt; // <===== NEW
}
}
これまた簡単です。でもまだ座標が更新されていないので正しくないので注意です。次は座標を更新していきましょう。
x を計算する
6.4 座標...、振り子の円弧の長さ やっと最後です!!ここまでくると次の座標
ただ、今回は 振り子の円弧の長さx について考えているので、
public class Pendulum : MonoBehaviour
{
// t+dt秒の計算をするためにt秒の情報が必要なので、取っておくための変数
float _x;
float _v;
void Start()
{
// 初期位置をいれておく。開始時の振り子の円弧の長さ。
_x = 1;
}
void FixedUpdate()
{
// [説明変数] 定数
float dt = Time.deltaTime;
float g = Physics.gravity.magnitude;
float l = 5;
// [説明変数] 現在の値
float x_t = _x;
float v_t = _v;
// 加速度(速度の時間微分)
float a = -g / l * x_t;
// dt秒先の速度を計算
_v = v_t + a * dt;
// dt秒先の座標を計算
_x = x_t + v_t * dt; // <===== NEW
}
}
また一行追加だけです。これで振り子の円弧の長さ
7. Unityでシミュレーションを可視化
7.1 Vector3のワールド座標を計算する
実際に球のGameObjectを動かすにはワールド座標
得られた振り子の円弧の長さ
いろいろ手段はありそうですが、運動方程式に
この
回転はクォータニオンを使います。Quaternion.AngleAxis()
を使うことで軸を指定して、好きな角度だけ回転するようなクォータニオンを手に入れられます。なので、振り子の支点を軸として
物理シミュレーションで得た 振り子の円弧の長さ
また、非常に申し訳ないですが、
必要であれば追記するか、別記事に補足として書きたいと思います。
Vector3 CalcBallPosition(float x, float l)
{
// 円弧の長さが分かっても座標が分からないので、運動方程式にあったθを求めて座標を計算する
// F = -m * g * sinΘ ≒ -m * g / l * x
// sinΘ = x / l
// Θ = arcsin(x / l)
float theta = Mathf.Asin(_x / l) * Mathf.Rad2Deg;
// クォータニオンを用いてθだけ初期位置から回転させる
var basePosition = l * Vector3.down;
var swingingPosition = Quaternion.AngleAxis(theta, transform.forward) * basePosition;
return swingingPosition;
}
7.2 可視化する
先ほど定義したメソッドで 振り子の円弧の長さ
デフォルトの球体のGameObjectの座標をこれで更新します。
また、振り子なので紐が欲しいですね。自分で動的Mesh生成してもよいですが、Unityには LineRenderer
というコンポーネントがあり、今回はこれで十分なので使わせてもらいます。
それらを加えた最終的なコードは以下になります。
using UnityEngine;
public class Pendulum : MonoBehaviour
{
[Header("Components")]
// 振り子の球
[SerializeField] Transform _ball;
// 振り子の紐
[SerializeField] LineRenderer _line;
[Header("Parameters")]
// 振り子の紐の長さ
[SerializeField] float _length = 5;
// t+dt秒の計算をするためにt秒の情報が必要なので、取っておくための変数
float _x;
float _v;
void Start()
{
// 初期位置をいれておく。開始時の振り子の円弧の長さ。
_x = 1;
}
void FixedUpdate()
{
// [説明変数] 定数
float dt = Time.deltaTime;
float g = Physics.gravity.magnitude;
float l = _length;
// [説明変数] 現在の値
float v_t = _v;
float x_t = _x;
// 加速度(速度の時間微分)
float a = -g / l * x_t;
// dt秒先の速度を計算
_v = v_t + a * dt;
// dt秒先の座標を計算
_x = x_t + v_t * dt;
var swingingPosition = CalcBallPosition(_x, l);
_ball.localPosition = swingingPosition;
var lineRootPosition = transform.position;
_line.SetPositions(new []{lineRootPosition, lineRootPosition + swingingPosition});
}
Vector3 CalcBallPosition(float x, float l)
{
// 円弧の長さが分かっても座標が分からないので、運動方程式にあったθを求めて座標を計算する
// F = -m * g * sinΘ ≒ -m * g / l * x
// sinΘ = x / l
// Θ = arcsin(x / l)
float theta = Mathf.Asin(_x / l) * Mathf.Rad2Deg;
// クォータニオンを用いてθだけ初期位置から回転させる
var basePosition = l * Vector3.down;
var swingingPosition = Quaternion.AngleAxis(theta, transform.forward) * basePosition;
return swingingPosition;
}
}
実際に動かしてみた様子
8. まとめ
- テイラー展開によって、ある
付近なら未来の値が予測できるa - テイラー展開によって、毎フレームで直近の座標を予測し、物理シミュレーションをした
今回はテイラー展開の
シンプレクティック法により、簡単な変更で精度が改善したりします。既に文章量が多いので割愛しますが、ほんの少し変更を加えるだけなのでもしかしたら追記するかもしれません。他にも改良オイラー陽解法だったり、他の方法などもあります。
おわりに
今回は振り子のシミュレーションをしました。減速するバネの運動方程式もあるので、衝突時の力から加速度を計算してCubeのScaleを動かしてみると、バウンドする気持ちいい物理的なアニメーションが作れるのではないかと思います。
終盤では、説明を省いてしまったものもいくつかありますが、本題の物理シミュレーション以外で話が膨らむと理解のカロリーが高くなるかと思い割愛させていただきました。
まだ自分自身も勉強中ですので、他の手法などについてもメモ代わりにまた学びながら書いていこうと思います。
Discussion