Mathematica では繰り返しのコードを書いたら負け」とは、以前の記事に矢吹さんがトラックバックで書かれていた言葉。まったく仰る通りなのだが、今回もまた繰り返しの話(つまり最初から負け...)。

何かの値を少しずつ変えながら同じ計算を繰り返し行なって様子を見たいことがある。このようなときに活躍するのが ver.6 で追加された Manipulate で、ちょっとした実験を行なうのが本当に容易になった。だが、Manipulate には向かないケースもある。1回の計算に数秒以上の時間がかかってしまう場合や、どこまで繰り返したら望んだ結果にたどり着くかが事前に分からないような場合だ。

簡単な例として、実数の逆行列を求める計算の速度が行列のサイズでどのように変わるかを見てみたいとする。行列サイズが大きくなれば計算時間が長くなるのはやる前から明らかだとしても、いったいどれくらいのサイズまで行列を大きくすれば、行列サイズ vs 計算時間のグラフが、見た目にそれらしい(自分が欲しい)ものになるかはやってみないと分からない(だいたい分かってるという突っ込みはなしでw)。

行列サイズを少し変化させるごとに行列サイズ vs 計算時間のグラフを新しく描画するというのは、ver.5.2 までのやり方。ループの回数だけグラフィックスがノートブックに追加されて、ノートブックファイルがすぐに何十メガバイトにもなってしまう。単に計算の進行度合いを見たかっただけなのに、どうしてこんな巨大ファイルを扱う羽目に...、というのは、古くからの Mathematica ユーザーの多くに経験のある感情ではないだろうか。ノートブックセルを扱う関数を駆使すれば、いちど描いたグラフのセルを差し替え描画することも可能ではあったが、少し面倒だった。

Mathematica 8 なら Dynamic を使うことで、いちど描いたグラフを動的に更新することができるので実に簡単だ。実数乱数の行列の逆行列算出にかかる時間を5回計測して平均値を出し、これを行列サイズを1つずつ大きくしながら行列サイズ vs 計算時間のグラフを更新していくというのは、次のように書ける。

data = {};
Dynamic[ListLinePlot[data, PlotRange -> All]]
For[n = 1, True, n++,
 times = Map[Function[{m}, First@Timing[Inverse[m]]],
   RandomReal[{-1, 1}, {5, n, n}]];
 AppendTo[data, {n, Mean[times]}]]

実行すると、グラフがみるみる更新されていく。この更新は、[評価] メニューの [評価を放棄] をするまで継続される。

この方法、つまり Dynamic[ListPlot[data]] を描いておいて、data を更新していく方法は、繰り返し計算で収束するような系の収束の状態をグラフで逐次確認したいときにも有効だ。もちろん、ListPlot でなくても良く、例えば昨年紹介した迷路作成は ArrayPlot で同じことをしている例だ。また、明示的なループを使わない FindFitNDSolve などでも、途中の評価点を StepMonitorEvaluationMonitor オプションで取り出すことができるので、これで data を更新すれば、どのように評価点を探索しているのかをリアルタイムで可視化して確認することが可能だ。

この、繰り返し計算中にリアルタイムで可視化を行なうというのは、規模の小さなシミュレーションにも使えるかもしれない。次の例は、粉体シミュレーションで使われる個別要素法をまねて、球と壁の接触箇所に垂直方向の弾性と粘性の力を算出し、ステップで粒子位置を更新するような計算をしたものだ。Graphics3DViewVertical オプションを Dynamic で指定すると、描画されたオブジェクトをマウスでぐりぐりするのに応じて、重力の方向が変わって球の転がりが変わる計算が無限に行なわれる。

dp = 0.5; rho = 0.9*^3; gravdir = {0, 0, 1};
pos = Module[{p, s, d, e, f, v, m = rho Pi dp^3/6, dt = 3.*^-2},
  p = v = {0, 0, 0};
  Function[{gd},
    s = Clip[p] - p; (* 壁との衝突 *)
    d = -0.5*^3 Sign[s]^2 v; (* 粘性力 *)
    e = 0.2*^6 Sign[s] Abs[s]^(3/2); (* 弾性力 *)
    f = d + e - m 9.8 Normalize[gd]; (* 力の合計 *)
    v += f/m dt; p += v dt]]; (* 速度を算出して座標更新 *)
Graphics3D[Sphere[Dynamic[pos[gravdir]], dp/2],
  PlotRange -> 1 + dp/2, SphericalRegion -> True,
  ViewVertical -> Dynamic[gravdir]]

この例では、各種パラメーターは見た目がそれっぽくなるように適当に決めたし、剪断方向の力や摩擦も組み入れていないので、これがシミュレーションとしてそのまま使えるものだと言うつもりは毛頭ない。だが、こんな風に計算をリアルタイムで可視化できるというのは、計算方法を検討する段階の小規模検算環境として、すごく便利だと思うのだけど、どうだろう。

※この記事の内容は執筆者の個人的見解で、ヒューリンクスによる公式情報ではありません。[免責事項]

トラックバック

この記事へのトラックバックURL
http://blog.hulinks.co.jp/cgi/mt/mt-tb.cgi/471
内容に対しての関連性がみられないものは削除する場合があります

コメント一覧

Dynamicの機能をうまく生かした面白い例題ですね。Mathematicaでは繰り返しを使わず関数型で書くのが常套ですが、Dynamicは動的に変化する対象を明示的に指定する必要があるので、計算の中間状態を明示的に持たない関数型とは相性悪い。でも、解はありそうな気がします。かなり昔から関数型の世界で動的な要素をどう扱うかという議論があり、古い論文などをひっくり返してみます。刺激的なブログ、感謝です。

skybird さん、いつもありがとうございます。何か新しい情報をご入手された折には、是非教えて下さい。できればかいつまんでw。

コメントの投稿

Emailアドレスは表示されません。は必須項目です。
ヒューリンクス取り扱い製品の内容や購入に関するお問い合わせはヒューリンクスサイト連絡先へお願いいたします。投稿前にその他の注意事項もご覧ください。

HULINKS サイトの新着情報