Mathematica で YouTube 動画を作る（Lorentz 方程式の数値解の可視化）

はじめに

吉野です．Mathematica でアニメーションを作るのが好きです．その一部は YouTube で公開しています．ここでは，以前に作った Lorentz 方程式の動画の作成手順を解説します． Mathematica で YouTube 動画を作ることがそれほど難しくないと思っていただければ幸いです．

Lorentz 方程式

Lorentz 方程式はカオスを示すことで有名な方程式です．

\begin{eqnarray*}
\frac{dx}{dt} &=&\sigma (y-x) \\
\frac{dy}{dt} &=& x(\rho-z)-y \\
\frac{dz}{dt} &=& xy-\beta z
\end{eqnarray*}

カオスの特徴のひとつに，時間が経過するにつれて初期値の僅かな違いがとても大きな違いになることが挙げられます．

Mathematica で数値解を求める

パラメータの値を固定して，初期値を指定すれば NDSolve 数値解を求めることができます．

lorentzEq = 
 With[{sigma = 10, rho = 28, 
   beta = 8/3}, {x'[t] == sigma (y[t] - x[t]), 
   y'[t] == x [t] (rho - z[t]) - y[t], 
   z'[t] == x[t] y[t] - beta z[t], x[0] == x0, y[0] == y0, 
   z[0] == z0}];
ansLorentz = 
 NDSolve[(lorentzEq /. {x0 -> 1, y0 -> 1, z0 -> 1}), {x, y, z}, {t, 
   0, 50}];
ParametricPlot3D[{x[t], y[t], z[t]} /. ansLorentz, {t, 0, 30}]

で，数値解をプロットすることができます．

初期値に敏感なことを示したい

「初期値に敏感であることを視覚的に示すには，微妙に異なる初期値についての数値解が十分に時間を経過すると異なる振る舞いをすることを示せばよい」と考えました．数学的にはリヤプノフ指数を使います．そこで，微妙に異なる初期値を20000ほど用意してその数値解を求め，相空間内での軌跡を表示するアニメーションを作成しました．

数値計算は以下のコードで行いました．初期値は， (0, 0, 1) を中心に半径 0.1 の球の中に入るように設定しました．コードは Table で 20000 回の繰り返しを行うだけの簡単なものです．

ansList = 
   Table[
    With[{r = RandomReal[{0, 0.1}], theta = RandomReal[{0, 2 Pi}], 
      phi = RandomReal[{0, Pi}]}, 
     NDSolve[(lorentzEq /. {x0 -> 1 + r*Cos[theta]*Sin[phi], 
         y0 -> 1 + r*Sin[theta]*Sin[phi], z0 -> 1 + r*Cos[phi]}), {x, 
       y, z}, {t, 0, 50}]
     ], {20000}]);

アニメーションは以下のコードで作成しました．表示に工夫したせいで，完成まではそれなりに時間がかかった記憶があります．でも，それほど難しいことはしていません．

grdata = Table[
   Graphics3D[
    Table[Sphere[{x[t], y[t], z[t]}, 0.2] /. ansList[[i]], {i, 
      20000}], PlotRange -> {{-30, 30}, {-30, 30}, {0, 60}}, 
    Boxed -> False, Background -> Black, Lighting -> "Neutral", 
    ImageSize -> 1200], {t, 0, 25, 0.02}];
Export["Lorentz20000-0-25-1240.mov", Take[grdata, 1240]];

結果を YouTube にアップロードしています（こちらからどうぞ）．20000 個の粒子は，はじめのうちは塊になって動いていますが，あるタイミングで細長くなり，最終的には別の動きをすることが確認できます．異なる軌道にあるはずの粒子がストレンジ・アトラクタの形を作るところなどはかなり面白いと思うのですが…．

是非とも Mathematica で数理アニメーションを作って公開してみてください．