6. グラフ作成とオイラー法の結果の可視化

赤木 ：前回、数字だけで、ちゃんと数値解がもっともらしい軌道になってい るかどうかみてなかったから、今回グラフ書いてみましょう。

学生C：いいですけど、何で書きます？ gnuplot とか？

赤木 ：それでもいいんだけど、プログラムの中から直接グラフ書けるほうが 便利よね。

学生C：うーん、まあ、そうですかね。データファイルからグラフでも、、、

赤木 ：まあ、今回、折角だから GR Framework ってのを使ってみましょう。これ、C で書いてあって Cからも Fortran からも使えてあと Python とか Julia からも 使えるっぽいので、 Crystal からも使って見れるかなと。

学生C： その辺よくわかってないですが、、、

赤木 ：まあ一緒にやってみましょう。GR のインストールはできたとして、 環境変数 GRDIR がインストールディレクトリをさしてるとすると、

 #include <stdio.h>
 #include <gr.h>
 
 int main(void) {
     double x[] = {0, 0.2, 0.4, 0.6, 0.8, 1.0};
     double y[] = {0.3, 0.5, 0.4, 0.2, 0.6, 0.7};
     //    gr_polyline(6, x, y);
     gr_axes(0.25, 0.25, 0, 0, 1, 1, -0.01);
     // Press any key to exit
     getc(stdin);
     return 0;
 }

 cc -I $GRDIR/include grsample.c -o grsample -L $GRDIR/lib -Wl,-rpath,$GRDIR/lib -lGR

 ./grsample

インストールは C ラ イブラリのインストールページから、Ubuntu なら gr-latest-Ubuntu-x86_64.tar.gz を落としてきて、それを どっか適当なディレクトリで tar xzf とかで展開するだけ。そうすると、grっ てディレクトリと、その下に libとか include とかのサブディレクトリができるの。そのディレクトリの名前 を GRDIR に設定して。bashなら

  export GRDIR=/foo/var/gr

  setenv  GRDIR /foo/var/gr

  export GKS_WSTYPE=x11

学生C：とりあえずできました。

赤木 ：そうねえ、今回とりあえず、作者氏が作った Crystal から GR の関数 を呼ぶためのライブラリ gr-crystalでごまかし てみて。

学生C：どうするんですかね？

赤木 ：そこに書いてあるから読んでやってみて。

学生C：えーと、 shard.yml になんか追加しろと。そこにある分だけでいいで すかね？

赤木 ：あ、最初になんかいるはず。以下でやってみて:

 name: shards
 version: 0.1.0

 dependencies:
    grlib:
     github: jmakino/gr-crystal
     branch: master

  > shards install
  Fetching https://github.com/jmakino/gr-crystal.git
  Installing grlib (0.1.0 at master)

赤木 ：じゃあ、そこの README.md にある通り、

   crystal lib/grlib/examples/grsample.cr

学生C： あ、なんかでました。

 require "grlib"
 GR.setwindow(0,1,0,1)
 GR.box
 GR.polyline([0.0, 0.2, 0.4, 0.6, 0.8, 1.0],  [0.3, 0.5, 0.4, 0.2, 0.6, 0.7])
 GR.setcharheight(0.05)
 GR.mathtex(0.3,0.2,"x^2+y^2")
 c=gets

  require "grlib"

その次からの GR. で始まるのが全部、 GR の機能を Crystal から 呼ぶためのライブラリ関数を使ってるの。関数の形や引数の意味は GR の Python API と同じだから、そっちみて？

学生C：ええと、、、 gr.setwindow ってのがありますが、これですか？

赤木 ：そう。

学生C： えーと、すみません、全然分からないんですが、、、

   setwindow establishes a window, or rectangular subspace, of world
   coordinates to be plotted. If you desire log scaling or
   mirror-imaging of axes, use the SETSCALE function.

赤木 ：あー、それ、 GR のベースになってる GKS っていう、コンピュータグ ラフィックスシステムの規格の用語ね。ISO標準にもなってるのよ。と、 そんなのはどうでもいいわね。この場合、2次元のグラフで、横軸と縦軸の 範囲を指定するの。この場合、 x軸の範囲もy軸の範囲も 0 から1まで、とい うことね。

次の box は、Cの gr_axes を呼んでるんだけど、ちょっと工夫してあって GRのデフォルトでは下と左にしか軸を書かないのを、上と右にも書くように しているの。あと、7個引数があって面倒くさいから、適当にデフォルト値 を設定して、それでよければそれで、としちゃってるの。

polyline は「折れ線」を引くの。x座標の配列とy座標の配列を引数にとるの ね。Cの gr_polyline は最初に点の数をいれたけど、Crystal の配列は 自分の要素の数を .size でもってるから、そっちにを使うようにしたのね。

学生C： x と y で要素数違ったらどうするんですか？

赤木 ： 小さいほう使うんじゃないかしら？やってみて。

setcharheight(0.05) は文字の高さを決めるもので、これはグラフを表示して いるウィンドウとかの全体が左下 (0,0)で右下が (1,1) になる座標系での高 さのはず。なので、0.05 とか小さな値なの。この座標系を規格化装置座標とか いったような気が。normalized device coordinates ね。

学生C：最初はどうなってるんですか？

赤木 ：GR のドキュメントにどっかに書いてあるから捜して。そういうの見つ けられるようになるのも大事だから。

学生C：はい、、、

赤木 ：最後に、 mathtex(0.3,0.2,"x^2+y^2") だけど、これは、 最初の2つが位置で、これもさっきの規格化座標のほう。あと、 3個目の文字列だけど、この関数だと LaTeX の数式モードの中に いるとしてコンパイルされて表示されるわ。

学生C：もっと普通なのはないんですか？

赤木 ：捜してみて。なんとか text みたいなのがあるはず。

学生C：はあ。

赤木 ：でね、 やってみて欲しいのは、前回の単振動の数値解で

• 解の軌跡を x-v平面に表示する
• 1周じゃなくて5周とか10周書いてみる
• ステップサイズもいくつかかえてみる

学生C： えーと、今日ですか？

赤木 ： それはおまかせ。 学生C： それでは、、、、

(しばらくあと)

作ってみました。プログラムは

 require "grlib"
 include Math
 include GR
 x=1.0; v=0.0; k=1.0
 n=ARGV[0].to_i
 norb=ARGV[1].to_i
 wsize=ARGV[2].to_f
 h = 2*PI/n
 p! h
 t=0
 setwindow(-wsize, wsize,-wsize, wsize)
 box
 setcharheight(0.05)
 mathtex(0.5, 0.06, "x")
 mathtex(0.06, 0.5, "v")
 
 (n*norb).times{
   xp=x
   vp=v
   dv = -x*k*h
   x+= v*h
   v+= dv
   t+= h
   GR.polyline([xp,x], [vp,v])
 }
 p! [x, v, t]
 p! [x-cos(t), v-sin(t)]
 gets

赤木 ：なかなかよい感じね。

学生C：これ、図にしてみると分かりますが円というからせんですね。一応、 ステップ数増やすと段々円に近くなりますが、、、

赤木 ：そうね。何故そうなるかわかる？本当は円になるはずよね？

学生C：えーと、そういう方法だから、とか、、、

赤木 ：それはそうなんだけど、君が何故そうなるかわかってる感じがあんま りなかったり。今の時刻の値を 、 次の時刻の値を、とすると、次の時刻の値はど うなるかしら？

学生C：えーと、

ではないかと。

赤木 ：そうね。で、これから、 i=0 の値からの一般解を計算できるでしょ？

学生C：？？？

赤木 ： をベクトル とすれば

はいいわね？

学生C：えーと、上の式と比べて、、、多分。

赤木 ：そうすると、ここにでてくる行列を A とすれば、

でしょ？なので、あとは を計算できればよくて、これは が の形で対角化できれば、 が になる、ってのは線型代数でやったわね？

学生C：そういわれるとそういう気もします。

赤木 ：まあだから、そういうふうに機械的に手を動かせば答がでる、っての が線型代数の偉大さなんだけど、ここはちょっと格好つけて複素数でやってみ て。 としたら は？

学生C：(n*10分ほど計算) だから、 です。

赤木 ：複素数を掛けることは回転+拡大・縮小と考えられる、ってのは聞いたことくらいはある?

学生C：あるかもしれませんが、、、

赤木 ：まあやってみれば、、、 を、 が実数として、 の形にしてみて。

学生C： ということですね。

でしょうか。

赤木 ：はい、よくできました。 の絶対値が で、 回転角は からでるわね。

学生C：はい、そうです。

赤木 ：じゃあ、初期値を とすれば は？

学生C： で、回転しなが ら絶対値が大きくなるということですね。

赤木 ：そう。 が小さくなれば、1ステップでの絶対値の増えかたは だから に比例して小さくなるけ ど、同じ時刻まで積分するのに必要なステップ数は に比例して増え るので、結局どれだけ本当の解からずれるかは に比例するわけね。

学生C：角度のほうもずれませんか？

赤木 ：もちろんずれるけど、じゃあそれは の何乗か、計算してみて。

学生C：じゃあこれは読者の皆様の練習問題ということで、、、、

赤木 ：あらあら。まあいいわ。これから実際に計算していこうという問題は もちろん、調和振動じゃなくてもっとややこしい非線型な問題なわけだけど、 数値計算法の振舞いは基本的に線型な問題に対して定義されてるの。

ここまでの話で、調和振動に前進オイラー法を適用すると、ある時刻 ま で積分した時の誤差が に比例する、と示せたわけじゃない？

学生C：えーと、ちゃんと証明になってましたっけ？

赤木 ：そこはちゃんと証明にしてみて。練習問題ね。

学生C：はあ、、、

6.1. 課題