宇宙の深淵からこんにちは、
\displaystyle\int {\rm d}\boldsymbol{r}\hat{\psi}^{\dagger}(\boldsymbol{r})ぽっぴんフレンズ\hat{\psi}(\boldsymbol{r})
です！今回は Unity で Lorentz 方程式の軌道を描画することを内容として記事にしたいと思います！(Lorentz 方程式について知らない方もいるかと思いますが、結果をチラ見せした方が早いでしょう。今回は下の画像のような曲線を描くことを目標にします！もしかしたら、どこかで見たことがあったりするのではないでしょうか？)
ということで、Lorentz 方程式及びアルゴリズムの紹介、そしてコードの実装を見ていきましょう！
Unityのバージョン：Unity 2018.4.6f1
こんな方にオススメ：
- Unity でゲーム制作以外のことをしてみたい
- 数値計算の結果をいい感じに視覚化したい
- ローレンツ方程式が大好き
Lorentz 方程式とは
Lorentz 方程式は、以下のように表されます。
\begin{align}
\frac{{\rm d}x}{{\rm d}t}&=-px+py\\
\frac{{\rm d}y}{{\rm d}t}&=-xz+rx-y\\
\frac{{\rm d}z}{{\rm d}t}&=xy-bz
\end{align}
$t$ は時間、$x, y, z$ は時間 $t$ の関数あり物体の座標、$p, r, b$ は実数の定数です。左辺は $x, y, z$ の時間微分、右辺は $x, y, z$ を変数とする関数となっています。
Lorentz 方程式は対流の研究のために流体力学の方程式を非常に簡単にしで出来たものらしいのですが、初期値鋭敏性という非常に興味深い特徴がみられる方程式として有名です。初期値鋭敏性とは、簡単に言うと「初めの状態(位置)のほんの僅かな誤差が後々大きな違いとなって現れる」ことです。物理現象を予測する上で初めの状態を測定することは非常に重要ですが、測定には当然誤差がついてきます。Lorentz 方程式ではその誤差の影響で後々の結果が大きく変わってきてしまう、すなわち現象の長期的な予測が予測が非常に困難となるのです！天気予報で遠い未来の天気を予測したいときほど予報が大きく外れるのと同じですね！
方程式といえば、普通は「解」がありますよね？連立方程式であれば片方の文字を消去するなどすれば $x$ と $y$ の値が求まります。 Lorentz 方程式にももちろん解があることが想像されますが、実はこの Lorentz 方程式、解を厳密に(手計算で)求めることができません！(な、なんだってー！！！) しかし、こののような場合のために、厳密な解ではなくそれっぽい解(数値解)を求める方法が存在します。詳しくは次節で取り扱いますが、この数値解を計算することで、手計算で解が出なくとも、方程式のおおまかな振る舞いを知ることができるのです！
#Runge-Kutta 法
これから Runge-Kutta 法というアルゴリズムについて紹介します。ぶっちゃけると、本節の内容を理解していなくても、次節で紹介するコードなどをパクってしまえば実装はできてしまいます。なので、この節は最悪読み飛ばしても構いません(本節も一生懸命書いたので、できれば読み飛ばさないでほしいです(迫真))。
Lorentz 方程式のような、関数の微分を含む方程式を微分方程式と呼びます。とくに Lorentz 方程式は時間 $t$ の関数が $x, y, z$ の3つあり、これらが互いに絡まりあっているので、連立微分方程式と呼ばれます。微分方程式は手計算で解を求められるものもありますが、そうでないものが圧倒的に多いのです。そして、前節でも述べたとおり、Lorentz 方程式も手計算では解を求められないものになります。そこで、コンピュータのパワーを借りることによて、「それっぽい解」を求める方法を用います。今回はいくつかある方法の中でも、精度が高い** Runge-Kutta 法**(厳密には4段4次 Runge-Kutta 法)というものを用います。
簡単のため、まずは時間 $t$ の関数が $x$ 1つのみである1階微分方程式( $x$ の $t$ による1階微分のみを含む方程式)を考えましょう。一般的な1階微分方程式は以下のように書かれます。
\frac{{\rm d}x}{{\rm d}t}=f(x, t)
左辺は $x$ の $t$ による1階微分、右辺は $x, t$ による何かしらの関数です。この方程式において、時刻 $t=0$ における $x$ の値 $x_0$ から、一定のアルゴリズムに従って次々と方程式の解を生成します。記事の趣旨上、詳しい説明は省き、アルゴリズムの内容のみを紹介します。
h\ は十分小さい正の実数、x_i\ は時刻\ t_i=ih\ における\ x\ の値とする。\\
時刻\ t_{i+1}=(i+1)h\ における\ x\ の値\ x_{i+1}\ を以下の手順で求める。\\
\begin{align}
\\
&k_1=h\times f(x_i, t)\\
&k_2=h\times f(x_i+\frac{k_1}{2}, t+\frac{h}{2})\\
&k_3=h\times f(x_i+\frac{k_2}{2}, t+\frac{h}{2})\\
&k_4=h\times f(x_i+k_3, t+h)\\
\\
&k=\frac{k_1+2k_2+2k_3+k_4}{6}\\
&x_{i+1}=x_i+k
\end{align}
以上が Runge-Kutta 法のアルゴリズムです。時刻 $t_0=0$ における $x$ の値 $x_0$ を自分で決めておけば、上記のアルゴリズムに従って、時刻 $t_1=h$ における $x$ の値 $x_1$, 時刻 $t_2=2h$ における $x$ の値 $x_2$, $\cdots$, 時刻 $t_i=ih$ における $x$ の値 $x_i$, $\cdots$ が次々に求まるという寸法です。この方法は精度が高く( $x(t+h)$ の $h$ についてのテイラー展開に対し $h^4$ の項まで一致する)、微分方程式の数値計算によく用いられます。
上記のアルゴリズムは $t$ の関数が $x$ の1つだけの場合のものでした。Lorentz 方程式は$t$ の関数が $x, y, z$ の3つあるので、少し改良を加えましょう。まず、$\vec{r}=(x, y, z)$ として $x, y, z$ を1つのベクトルにまとめてしまいましょう。$\vec{r}$ は言うなれば $x, y, z$ を座標とする位置ベクトルになります。そうすると、$x, y, z$による連立1階微分方程式は以下のように表されます。
\frac{{\rm d}\vec{r}}{{\rm d}t}=\vec{f}(\vec{r}; t)\\
\begin{align}
ただし、&\vec{r}=(x, y, z),\\
&\vec{f}(\vec{r}; t)=(f_x(x, y, z; t), f_y(x, y, z; t), f_z(x, y, z; t))
\end{align}
左辺は $\vec{r}$ の各成分を微分したもの、右辺は $\vec{r}, t$ を変数とするベクトル関数です。ベクトル関数 $\vec{f}(\vec{r}; t)$ は、$x, y, z$ 成分それぞれが異なる形の関数 $f_x, f_y, f_z$ となっています。この微分方程式の Runge-Kutta 法は以下のようになります。
h\ は十分小さい正の実数、\vec{r}_i\ は時刻\ t_i=ih\ における位置とする。\\
時刻\ t_{i+1}=(i+1)h\ における位置\ \vec{r}_{i+1}\ を以下の手順で求める。\\
\begin{align}
\\
&\vec{k}_1=h\times \vec{f}(\vec{r}_i, t)\\
&\vec{k}_2=h\times \vec{f}(\vec{r}_i+\frac{\vec{k}_1}{2}, t+\frac{h}{2})\\
&\vec{k}_3=h\times \vec{f}(\vec{r}_i+\frac{\vec{k}_2}{2}, t+\frac{h}{2})\\
&\vec{k}_4=h\times \vec{f}(\vec{r}_i+\vec{k}_3, t+h)\\
\\
&\vec{k}=\frac{\vec{k}_1+2\vec{k}_2+2\vec{k}_3+\vec{k}_4}{6}\\
&\vec{r}_{i+1}=\vec{r}_i+\vec{k}
\end{align}
以上のアルゴリズムによって、時刻 $t_0=0$ における位置 $\vec{r}_0$ をはじめに決めておけば、方程式の解となる点が次々と求まります。
以上が、今回用いる Runge-Kutta 法の説明でした。少々難しかったでしょうかね(汗)。次節からいよいよUnity上で実装していきます！
#実装
いよいよ、Unity上で Lorentz 方程式を実装していきます！前節の Runge-Kutta 法が良くわからなくても、以下の内容に沿っていけば誰でも実装できますので、頑張っていきましょう！
##Unity上での準備
はじめに、Unity で新しいプロジェクト(3D)を開きます(Unityがない人は諦めるか今すぐインストールしてください！)。HierarchyのCreateからCreate Emptyを押し、空の GameObject を作ります。名前は Point としましょう。また、InspectorにあるTransformで Object の位置を設定します(これが初めの位置になります)。今回は $(x, y, z)=(0, 10, 1)$ としましょう。この GameObject に色々なものを付け足すことで、Lorentz 方程式に従って動きながら、軌道が線として残るようにします。
カメラの位置と向きも調整しておきましょう。軌道はかなり大きくなるので、カメラも割合遠くに置きます。HierarchyのMain Cameraを押し、InspectorにあるTransformで位置を調整します。今回は $(x, y, z)=(0, 25, -60)$ としましょう。
これで最低限の準備はできました！次はいよいよスクリプトを書いていきます！
Lorentz 方程式の実装
いよいよ Lorentz 方程式を Unity の C# スクリプトで実装していきますが、ここで一つ面倒な問題があります。それは、Unity での座標軸の取り方と数学・物理学での座標軸の取り方が異なるということです。数学・物理では、水平面に $x-y$ 平面が置かれ、鉛直方向に $z$ 軸が置かれますが、Unity では水平面に $x-z$ 平面が置かれ、鉛直方向に $y$ 軸が置かれます。すなわち、$y$ 軸と $z$ 軸が入れ替わった感じになります。最初に紹介した Lorentz 方程式の座標は数学・物理学での慣習に則っているので、Unity 上で実装する際に元の式の $y$ と $z$ を入れ替えて式を少し変えます(入れ替えなくても動くには動くのですが、軌道が横倒しになった形になります)。以下が、Unity 用に $y$ と $z$ を入れ替えた式になります。
\begin{align}
\frac{{\rm d}x}{{\rm d}t}&=-px+pz\\
\frac{{\rm d}y}{{\rm d}t}&=xz-by\\
\frac{{\rm d}z}{{\rm d}t}&=-xy+rx-z
\end{align}
右辺だけでなく左辺の微分のところも入れ替えることに注意してください。これで準備は整いました。
さて、いよいよコードを書いていきましょう！Unity のProjectでCreateを押し、C# Scriptを生成します。スクリプト名は LorentzEquation としましょう。スクリプトを開いて以下のコードを書いていきます。
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public class LorentzEquation : MonoBehaviour
{
//Lorentz 方程式内の定数
float p = 10.0f;
float r = 28.0f;
float b = 8.0f / 3.0f;
//Runge-Kutta 法で用いる、時間の刻み幅
float h = 5e-3f;
//毎フレームごとに Runge-Kutta 法で GameObject の位置をアップデート
void Update()
{
RungeKutta();
}
//Runge-Kutta法の関数
void RungeKutta()
{
Vector3 k1, k2, k3, k4;
Vector3 position, deltaPosition;
position = this.transform.position;
k1 = h * Lorentz(position);
k2 = h * Lorentz(position + k1 / 2);
k3 = h * Lorentz(position + k2 / 2);
k4 = h * Lorentz(position + k3);
deltaPosition = (k1 + 2 * k2 + 2 * k3 + k4) / 6;
this.transform.position += deltaPosition;
}
//Lorents 方程式の関数
Vector3 Lorentz(Vector3 position)
{
float x = position.x;
float y = position.y;
float z = position.z;
// dx/dt=-px+pz
// dy/dt=xz-by
// dz/dt=-xy+rx-z
float xdot = p * (-x + z);
float ydot = x * z - b * y;
float zdot = -x * y + r * x - z;
return new Vector3(xdot, ydot, zdot);
}
}
}
// Lorentz 方程式の関数
Vector3 Lorentz(Vector3 position)
{
float x, y, z;
float xdot, ydot, zdot;
x = position.x;
y = position.y;
z = position.z;
// dx/dt=-px+pz
// dy/dt=xz-by
// dz/dt=-xy+rx-z
xdot = p * (-x + z);
ydot = x * z - b * y;
zdot = -x * y + r * x - z;
return new Vector3(xdot, ydot, zdot);
}
}
このコードを先ほど生成した GameObject にひも付すれば、Object がローレンツ方程式に従って動いてくれるようになります。
コード内の関数の説明をしていきます。まずは Lorentz 方程式の関数から行きましょう。方程式内の定数はコードの冒頭で $p=10, r=28, b=8/3$ に設定します。
//Lorentz 方程式内の定数
float p = 10.0f;
float r = 28.0f;
float b = 8.0f / 3.0f;
GameObject の位置を引数としてローレンツ方程式の右辺を計算し、方程式の左辺を返す関数を作ります。はじめに Vector3 型の引数の $x, y, z$ 成分を x, y, z にそれぞれ格納し、次に方程式の右辺の計算結果を xdot, ydot, zdot にそれぞれ格納します。最後に xdot, ydot, zdot を成分とする Vector3 型の変数を返り値として返します。これで Lorentz 方程式の関数の完成です。
// Lorentz 方程式の関数
Vector3 Lorentz(Vector3 position)
{
float x, y, z;
float xdot, ydot, zdot;
x = position.x;
y = position.y;
z = position.z;
// dx/dt=-px+pz
// dy/dt=xz-by
// dz/dt=-xy+rx-z
xdot = p * (-x + z);
ydot = x * z - b * y;
zdot = -x * y + r * x - z;
return new Vector3(xdot, ydot, zdot);
}
次に、Runge-Kutta 法で GameObject の位置をアップデートする関数について説明します。時間 $t$ の刻み幅 $h$ はコードの冒頭で設定しております。$h$ を大きくとると線が荒くなってしまいますが、小さくとりすぎると GameObject の動きが遅すぎになってしまいます。今回は $h=0.005(=5e-3)$ としました。(※ $h$ として1フレームの時間 Time.deltaTime を用いたい場合は関数の中で $h$ を設定する必要があります。たぶん。)
//Runge-Kutta 法で用いる、時間の刻み幅
float h = 5e-3f;
Runge-Kutta 法の関数は、毎フレームごとに GameObject の位置を取得し、Runge-Kutta 法のアルゴリズムに従って位置の変化量を計算し、それを位置に加算することで GameObject の位置を更新します。まず、Runge-Kutta 法で用いる Vector3 型変数 k1, k2, k3, k4 と、そのフレームでの位置を格納する変数 position 及び Runge-Kutta 法の計算結果を格納する変数 deltaPosition を宣言します。次に、position に GameObject の位置を格納します。そして Runge-Kutta 法のアルゴリズムに従って k1, k2, k3, k4, 及び deltaPosition を計算し、最後に GameObject の位置に deltaPosition を加算します。
//Runge-Kutta法の関数
void RungeKutta()
{
Vector3 k1, k2, k3, k4;
Vector3 position, deltaPosition;
position = this.transform.position;
k1 = h * Lorentz(position);
k2 = h * Lorentz(position + k1 / 2);
k3 = h * Lorentz(position + k2 / 2);
k4 = h * Lorentz(position + k3);
deltaPosition = (k1 + 2 * k2 + 2 * k3 + k4) / 6;
this.transform.position += deltaPosition;
}
上の関数は1フレーム分の処理なので、これを void Update() 関数内に書くことで毎フレーム処理を行ってくれるようになります。
//毎フレームごとに Runge-Kutta 法で GameObject の位置をアップデート
void Update()
{
RungeKutta();
}
これでスクリプトの実装ができました！念のため、本節冒頭で作った GameObject にスクリプトをひも付することを忘れないでください。次は GameObject の軌道を描く方法について紹介していきます！
##軌道の描画の実装
前小節で、GameObjectを Lorentz 方程式に従て動かすためのスクリプトを作りましたが、これだけでは Object の軌道が見えません。ここでは、Object の軌道を表示させるための設定をしていきましょう！
Unity には線を描くためのコンポーネントとして Line Renderer というものがあります。まずはこれを最初に作った GameObject に追加しましょう。Hierarchy に先ほど生成した Point があるのでそれを押し、InspectorにあるAdd Componentを押します。検索窓が出てくるので、そこで「Line Renderer」と検索してコンポーネントを追加します。コンポーネントを追加したら、線の太さを変えていきます。widthという名前で謎のグラフがある個所があります。左上のボックスの数字を変えると、線の太さを調整する上限を変えられます。デフォルトで1.0になっているはずですが、これだと太すぎるので、今回は数字を0.15にし、赤い線をマウスで一番上まで引っ張ります。これで線の太さが0.15になります。
次に、線を描くスクリプトを作っていきます。ProjectからC# Scriptを生成し、以下のコードを書きます。
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public class DrawLine : MonoBehaviour
{
//オObject の軌跡を描画する
LineRenderer line; // LineRenderer Componentを受ける変数
int count; // 線の頂点の数
// Start is called before the first frame update
void Start()
{
line = GetComponent<LineRenderer>(); // LineRenderer Componentを取得
}
void Update()
{
count += 1; // 頂点数を１つ増やす
line.positionCount = count; // 頂点数の更新
line.SetPosition(count - 1, this.transform.position); // Object の位置情報をセット
}
}
1フレームごとに線を描くための頂点を増やしていき、増やした頂点に GameObject の位置をセットすることで線を描くことができます。尚、頂点の番号は 0, 1, 2, ... のようにゼロから振られるので、新しく追加される点の番号は $点の総数-1$ となることに注意してください(たとえば100個の点で描かれる曲線の100個目の点の番号は99番)。
ここまでで、既に Lorentz 方程式の軌道を描画する準備は整いました！では、いざ実行してみましょう！
#実装結果
準備が整ったので、実行してみましょう！...とその前に、以下の点を今一度確認してみてください。
- GameObject に LorentzEquation のスクリプト、DrawLine のスクリプト及び LineRenderer がひも付されているか。
- LineRenderer の線の太さは適切か(オススメ：0.15)。
- LorentzEquation のスクリプト、DrawLine のスクリプトの内容に間違いはないか。
- GameObject の位置は適切か(オススメ：$(0, 10, 1)$ )。
- カメラの位置は適切か(オススメ：$(0, 25, -60)$ )。
何かしらのエラーがあればConsoleに表示されるので、エラーに従って直してください。
すべて確認できたら、いざ実行してみましょう！以下のように動くはずです。(よかったらTwitterのフォローお願いします！)
UnityでLorentz方程式の軌道を描画してみました。#Unity#ローレンツアトラクタ pic.twitter.com/SxJC76hU5l
— ∫drψ†ぽっぴんフレンズψ (@dr30793443) October 11, 2019
しばらくすると冒頭で載せた画像のようになります。
圧巻ですね。。。
#おまけ：カメラの移動
さらなる圧巻を求めて、軌道のまわりを回転しながら眺められるように、カメラに動きを付け足していきましょう！ProjectでC# Scriptを生成し、名前を CameraScript としましょう。そして、以下のコードを書いていきます。
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public class CameraScript : MonoBehaviour
{
public float omega = 0.30f; //半径
public float A = 60.0f; //角速度
public float height = 25.0f; //高さ
float theta; //回転角
float sinTheta; //回転角のsin
float cosTheta; //回転角のcos
void Update()
{
theta = omega * Time.time; //角度=時刻*角速度
sinTheta = Mathf.Sin(theta); //回転角のsinを計算
cosTheta = Mathf.Cos(theta); //回転角のcosを計算
//カメラの位置の代入
this.transform.position = new Vector3(A * sinTheta, height, -A * cosTheta);
//カメラの方向の代入
this.transform.forward = new Vector3(-sinTheta, 0.0f, cosTheta);
}
}
はじめに、回転の半径(A)、角速度(Omega)、および回転する位置の高さ(height)を決めます。変数の前にpublicを付けることで、カメラのInspectorから変数をいじれるようになりますが、付けなくても動きはします。
public float omega = 0.30f; //半径
public float A = 60.0f; //角速度
public float height = 25.0f; //高さ
各フレームにおけるカメラの回転角とその $\sin, \cos$ を格納する変数も用意しておきます。
float theta; //回転角
float sinTheta; //回転角のsin
float cosTheta; //回転角のcos
そして、以下の部分でカメラの位置と向きを毎フレーム更新していきます。はじめに theta, sinTheta, cosTheta の値を計算します。最初に計算してしまうことで、$\sin, \cos$ を使いたいたびに Mathf.sin 関数で計算をする手間が省けます。時刻 $t$ におけるカメラの位置は $r = (A\sin{\omega t}, h, -A\cos(\omega t))$ にします。こうすることで、$(x, y, z)=(0, 25, -60)$ からスタートし、上から見て反時計回りに回転してくれます。カメラの向きは this.transform.forward というベクトルをいじることで調整でき、$(-\sin{\omega t}, 0, \cos(\omega t))$ とすることで常に円の中心を向いてくれます。
void Update()
{
theta = omega * Time.time; //角度=時刻*角速度
sinTheta = Mathf.Sin(theta); //回転角のsinを計算
cosTheta = Mathf.Cos(theta); //回転角のcosを計算
//カメラの位置の代入
this.transform.position = new Vector3(A * sinTheta, height, -A * cosTheta);
//カメラの方向の代入
this.transform.forward = new Vector3(-sinTheta, 0.0f, cosTheta);
}
以上のコードによって、カメラが高さ height=25 にある半径 A=60 の円周上を、角速度 omega=0.30 で、円の中心を向きながら回転してくれます。以下がその結果です。
カメラを回転させた場合です#Unity#ローレンツアトラクタ#プログラミング初心者 pic.twitter.com/efsVBDl162
— ∫drψ†ぽっぴんフレンズψ (@dr30793443) October 11, 2019
より臨場感が増しましたね！
#まとめ
今回は Lorentz 方程式を Unity で実装し、軌道を描画してみました。皆様、いかがでしたか？個人的には思ったより簡単にできたなという印象だったので、Unity が使える方は是非試してみてください！また、この記事の続きとして、複数の点を少しだけずらして出発させ、初期値鋭敏性が現れる様子を見るというのもやる予定なので、楽しみにしておいてください！
もしこの記事が面白いと感じていただけたのならば、是非いいねやコメントをお願いします！また、式や文章の間違い、記事の書き方についてのアドバイス、その他取り上げてほしい題材のリクエスト等あれば遠慮なくコメントしてください！それでは、また次の記事でお会いしましょう！