5. 運動方程式1 (2020/1/24)

5.1. 背景説明的導入

赤木 ：じゃあ、まずは、コンパイラとかのインストール 1 からクラスとかの解説 4 までやってみて。で、わかんないとことかあったら教えて。

学生C：じゃあちょっとやってみます。

5.2. 一次元調和振動子

学生C：4 までやってみました。

赤木 ：どうだった？難しい？

学生C：難しくないといえば嘘になりますというか、、、なんというか、だからなに？みたいな、、、

赤木 ：そうよね。なので、この辺からもうちょっと楽しいことをしましょうよ、というわけ。

学生C：楽しいんですか？

赤木 ：まあ人によっては。

学生C：うーん、、、

赤木 ：まあしばらく付き合って。というわけで、これやってみて:

一次元調和振動子

1 初期条件

2 それぞれの方法で、 での誤差を、時間刻みの関数として グラフ化し、時間の何次になっているか、何故そうなるか議論せよ。

なお、ここでは、数値計算パッケージとかソフトウェアを使うのではなく、自 分でプログラムを書くこと。

赤木 ：前進オイラー法とかはそれぞれどういうのか知ってる？

学生C：多分、、、まず前進オイラー法ですね。常微分方程式の初期値問題で すから、方程式が

赤木 ：うーんと、、、多分。2次ルンゲクッタは？

学生C：同じ微分方程式として、

赤木 ：まあプログラム書いてみて上手くいけばあってるはずよね。前進オイラー法の原理は？

学生C：教科書的には、もとの常微分方程式の解が存在していてそのテイラー 展開も存在しているとすると、

赤木 ：じゃあ、それで計算したものが、この問題の、 まで計算した時に近似解になるのはどうして?

学生C：え、どうしてって、その、数値解だから近似解なのでは？

赤木 ：うーん、、、そうじゃなくて、、、じゃあまずはオイラー法だけでいいから実際にプログラム作って計算してみて？

学生C：はあ、、、やってみます。

(一週間後)

学生C：こんな感じです。

 #!/usr/bin/env crystal
 include Math
 x=1.0; v=0.0; k=1.0
 h=ARGV[0].to_f*PI
 p! h
 t=0
 while t< PI*2 - h/2
   dv = -x*k*h
   x+= v*h
   v+= dv
   t+= h
 end
 p! [x, v, t]
 p! [x-cos(t), v-sin(t)]

 gravity> ./harmonic-euler.cr 1e-3
 h # => 0.0031415926535897933
 [x, v, t] # => [1.0099184201712226, 2.0875749689053084e-5, 6.283185307179335]
 [x - (cos(t)), v - (sin(t))] # => [0.009918420171222575, 2.0875749940652505e-5]
 gravity> ./harmonic-euler.cr 1e-4
 h # => 0.0003141592653589793
 [x, v, t] # => [1.0009874475970644, 2.069126097899751e-7, 6.283185307177473]
 [x - (cos(t)), v - (sin(t))] # => [0.0009874475970643726, 2.0691472301136493e-7]
 gravity> ./harmonic-euler.cr 1e-5
 h # => 3.1415926535897935e-5
 [x, v, t] # => [1.000098700914577, 2.0672973102662693e-9, 6.283185307165095]
 [x - (cos(t)), v - (sin(t))] # => [9.870091457697683e-5, 2.0817890742914563e-9]
 gravity> ./harmonic-euler.cr 1e-6
 h # => 3.141592653589793e-6
 [x, v, t] # => [1.0000098696530915, 2.062700127729346e-11, 6.283185307154872]
 [x - (cos(t)), v - (sin(t))] # => [9.86965309146548e-6, 4.5341698913228974e-11]

学生C：では最初から。最初の

  #!/usr/bin/env crystal

赤木 ：うーん、そうねえ、でも、繰り返して実行するなら、 スクリプトじゃ ないからちゃんとコンパイルしたほうがいいわ。

  crystal build harmonic-euler.cr

  crystal build --release harmonic-euler.cr

学生C：あ、そういえばそうでした。では最初の行は今後なしで。

赤木 ：はい。次の行は？

学生C：

 include Math

赤木 ：Math は言語というか文法的にはどういうもの？

学生C：module ですね。class と module がどう違うのかよくわかってないで すが、module の中で定義したメソッドは、 class の中で include すると そのクラスのメソッドになるし、 class の中じゃなくて外だと普通の関数に なるようです。

赤木 ：そうね。次は x, v, k の値ね。それはいいとしてその次は？

学生C： ARGV は、コマンドラインで与えた引数が順番に入る配列です。これ は C言語の argv ですが、0から引数なのは Ruby と同じです。 PI は モジュール Math で定義されている円周率で、 2まで積分する ので入力した数値に円周率かけたものを実際の h にしました。

赤木 ：了解。あと残りは？

学生C：数値積分の本体は

  while t< PI*2
    dv = -x*k*h
    x+= v*h
    v+= dv
    t+= h
  end

    x+= v*h
    v-= x*k*h

時刻も0から加算していって、 2*PIまできたら止めるとしました。あとは

 p! [x, v, t]
 p! [x-cos(t), v-sin(t)]

最後に実行結果ですが、hを10倍づつ変えて、誤差をだすと大体1桁づつ小さく なるので、ちゃんとできてるんではないかと思いますが、どうでしょう？

赤木 ：えーと、そうね、一見いいように見えるんだけど、、、t の値は大丈 夫？

学生C：大丈夫というと？

赤木 ：ちゃんと になってる？

学生C：え、なってるでしょ？値が 6.2863268998329245、、、あれ？ が 3.141592 として 2倍は 6.283184、、、あれれ？4桁めから違いま すね。ちゃんと計算すると、、、折角だから電卓に crystal 使ってみます。

 gravity> crystal eval "p 6.2863268998329245/(Math::PI)/2"
 1.00049999999996

 gravity> crystal eval "p 1.0/(6.2863268998329245/(Math::PI)/2-1)"
 2000.0000001600924

赤木 ：そう。だから、1ステップ余計に計算しちゃってるの。どうしてかわかる？

学生C：本当は2000ステップちょうどで終わらないといけないんですが、2000 ステップでは にならないということですね。浮動小数点の足し算で 誤差があるから、、、

赤木 ：そう。どう直すのがいい？

学生C： 比較する値をちょっと小さいめに、例えば h/2 をひいておくとか、、、

赤木 ： それでもいいわね。もうひとつは、はじめから、何ステップで積分す るかのほうをいれて、hのほうをあとで決めるのね。そっちでやってみて。

学生C： times を使ってみました。n.times{色々} で 「色々」を n 回繰り返します。

 include Math
 x=1.0; v=0.0; k=1.0
 n=ARGV[0].to_i
 h = 2*PI/n
 p! h
 t=0
 n.times{
   dv = -x*k*h
   x+= v*h
   v+= dv
   t+= h
 }
 p! [x, v, t]
 p! [x-cos(t), v-sin(t)]

 gravity> crystal build harmonic-euler-v3.cr
 gravity> ./harmonic-euler-v3 1000
 h # => 0.006283185307179587
 [x, v, t] # => [1.0199349143076457, 8.432969374211775e-5, 6.283185307179473]
 [x - (cos(t)), v - (sin(t))] # => [0.019934914307645712, 8.432969385516134e-5]
 gravity> ./harmonic-euler-v3 10000
 h # => 0.0006283185307179586
 [x, v, t] # => [1.0019758699537742, 8.284675641517022e-7, 6.28318530718043]
 [x - (cos(t)), v - (sin(t))] # => [0.0019758699537741897, 8.284667206271328e-7]
 gravity> ./harmonic-euler-v3 100000
 h # => 6.283185307179586e-5
 [x, v, t] # => [1.0001974115707355, 8.269974924947398e-9, 6.283185307175863]
 [x - (cos(t)), v - (sin(t))] # => [0.0001974115707354951, 8.273698413812141e-9]
 gravity> ./harmonic-euler-v3 1000000
 h # => 6.283185307179587e-6
 [x, v, t] # => [1.0000197394036099, 8.271237919989742e-11, 6.283185307155417]
 [x - (cos(t)), v - (sin(t))] # => [1.9739403609886352e-5, 1.0688173528613706e-10]

5.3. 課題