1 調和振動

学生B：赤木さんこんにちは。お世話になります。

赤木： あら、こちらこそ、引き受けてくれてありがとう。

学生B：あ、いえ。これからシミュレーションのプログラムをなんとかしないといけないみたいなので、ちょうどいい機会だと思います。

赤木： そう。どこから始めるのがいいかしら？学部の講義で数値計算とかやってるわよね？ルンゲ・クッタとか線型多段階法ってわかる？

学生B：ええと、すみません、どちらも聞いたことはあるんですが、、、

赤木： うーん、しょうがないわね。オイラー法は？

学生B：あ、それはおぼえてます。

赤木： ちょっとどういうものか説明してみて。

学生B：ええと、確か、

$$\frac{dx}{dt} = f(x,t)$$ (1)

みたいな常微分方程式があったとして、

$$x(t+\Delta t) = x(t) + \Delta t f(x(t),t)$$ (2)

というふうに、 時刻 $t$ での解から時刻 $t+\Delta t$ の解を求める方法で す。

赤木： なんだか記号の使い方が豪快だけど、基本的にはそういうことね。 図1みたいに、時刻 $t_i$ での近似解を $x_i$ とすれば、次 の時刻 $t_{i+1} = t_i + \Delta t$ での近似解が

$$x_{i+1} = x_i + \Delta t f(x_i,t_i)$$ (3)

で求められるということね。

Figure 1: オイラー法の概念

図は $x$が1次元の場合だけど、多変数になっても原理は同じね。とにかく今 のところでの真の解の接線を引くわけ。

でも、あれよね、図に細い線で書いてあるのが「本当の」解 のつもりなわけね。そうすると、もちろんオイラー法で求まった解は本当の解 とは違うわけじゃない。どれくらい違うのかしら？数学としてちゃんとやるに は収束性とかが問題になるけど、これ割に面倒だからとりあえず簡単な例でプ ログラム書いて調べてみて。ええと、こうしましょう。

問題:以下の連立微分方程式

$$\frac{dx}{dt} = v$$ (4)

$$\frac{dv}{dt} = -kx$$ (5)

を、初期条件 $x_0 = 1, v_0 = 0$ から 1 周期オイラー法で数値解を求めて、1周期後の誤差を時間刻みの関数として求めてみよ。

学生B：時間刻みはどれくらいの範囲でやるんですか

赤木： そういうのは自分で考えてみてね。

学生B：(なんか適当な人だ、、、)はあ、、、

赤木： じゃあ、よろしく。結果はグラフに整理してね。

学生B：わかりました、、、じゃあ、出来たらまたきますね。

しばらくたっても Bが戻ってくる気配がないので、探しにいく赤木。

赤木： ねえ、まだ？

学生B：あ、、、その、プログラム書いて答もでることはでるんですが、どう もおかしいんで、直そうとしてるんですが、、、

赤木： どうおかしいの？

学生B：ええと、精度が良過ぎるような気がするんです。ちょっと実行結果を おみせします。

<yebisu:/home/makino/papers/Interface>euler
Enter k, nsteps :1 10
k, period, nsteps, dt = 1 6.28319 10 0.628319
0.628318530717958623 -0.041597661808345579 -0.113783383682847572
<yebisu:/home/makino/papers/Interface>!!
euler
Enter k, nsteps :1 100
k, period, nsteps, dt = 1 6.28319 100 0.0628319
0.0628318530717958679 -3.30347425558861194e-05 -0.00103451230852307213
<yebisu:/home/makino/papers/Interface>!!
euler
Enter k, nsteps :1 1000
k, period, nsteps, dt = 1 6.28319 1000 0.00628319
0.00628318530717958661 -3.25234124386782758e-08 -1.03355224671600188e-05
<yebisu:/home/makino/papers/Interface>

学生B：最初にいれてるのは $k$ で、これはいつでも 1 です。で、次にいれ てるのが1周期を何ステップでやるかで、これを 10, 100, 1000 と10倍ずつ変 えてみたんですが、、、 $x$ の誤差は 3 桁づつ減るし、$v$ の誤差も2桁づ つ減ってるんです。これって変じゃないですか？でも、プログラムいくらみてもおかしいところはないみたいだし、、、

赤木： まあ、とりあえずプログラム見せてくれる？

学生B：ええと、これです。

pt

//
// Euler.C
//
// J. Makino 2002/06/02
#include  <stdlib.h>
#include  <math.h>
#include  <iostream>

using namespace std;

int main()
{
    double x,v,t;
    double k;
    int nsteps;
    cerr << "Enter k, nsteps :";
    cin >> k >> nsteps;
    double period = sqrt(1.0/k)*M_PI*2;
    double dt = period/nsteps;
    cerr << "k, period, nsteps, dt = " << k << " "
         << period << " "
         << nsteps << " "
         << dt << endl;

    x=1;
    v=0;
    t=0;
    for(int i = 0;i<nsteps;i++){
        v += -dt*k*x;
        x += dt*v;
        t += dt;
        cout.precision(18);
        cout << t << " " << x << " " << v << endl;
    }
    cout << dt << " " << x-1 << " " << v  << endl;
}

赤木： ちょっと見せて。(数十秒後)ああ、なるほどね。それはそうなるわね。

学生B：ええっ？なんかおかしいことしてます？

赤木： (こいつもなかなか天然ではないか。。。)そうねえ、あら、ステップ毎に結果を書き出すこともできるのね。まず、最初のステップがちゃんとできてるかどうかってみた？

学生B：え、1周期でちゃんと戻ってきてるし、グラフにも書いてみてコサインの形になっているのは確認しましたけど。

赤木： そうじゃなくて、数字自体がちゃんとでてるかどうかよ。オイラー 法の式なんか簡単なんだから、電卓でも bcでも使って計算できるでしょ？

学生B：えーと、じゃあ、 $\Delta t = \pi/5$ の時でやってみます。最初は $x_0 = 1, v_0 = 0$ ですから、ええと、 $x_1 = 1, v_1 = -\pi/5 = .62831853$ くらいですね。

赤木： 「くらい」ってのはなによ？まあいいわ、で、プログラムの答は？

学生B：ちょっと待って下さい、、、ええと、

euler
Enter k, nsteps :1 10
k, period, nsteps, dt = 1 6.28319 10 0.628319
0.628318530717958623 0.605215823956425703 -0.628318530717958623
......

あれ、$x$が1じゃない、、、え、どうしてですか？

赤木： どうしてって、あなたのプログラムがそうなってるからそういう答になるんじゃないの？プログラムではどう計算してるの？

学生B：ええと、まず $v$を計算して、それから $x$を、、、あ、そういうこ とか、 新しい $x$ を計算するのに、$v$ の新しい値を使っているのが違うわ けですね。 とりあえず直して、、、

	double dv = -dt*k*x;
	double dx = dt*v;
	v += dv;
	x += dx;

こんなのでどうでしょう。

赤木： まあ、よさそうよね。動かしてみたら？

学生B：

<yebisu:/home/makino/papers/Interface>euler1
Enter k, nsteps :1 10
k, period, nsteps, dt = 1 6.28319 10 0.628319
0.628318530717958623 1 -0.628318530717958623
......

と、最初は今度はいいみたいです。

赤木： じゃあ、1周期やってみて。

学生B：では、ステップ毎に書くのをやめて、、、

<delaware:/usr2/makino/papers/Interface>!eu
euler1
Enter k, nsteps :1 10
k, period, nsteps, dt = 1 6.28319 10 0.628319
0.628318530717958623 3.12658511680165985 3.29196073410105283
<delaware:/usr2/makino/papers/Interface>euler1
Enter k, nsteps :1 100
k, period, nsteps, dt = 1 6.28319 100 0.0628319
0.0628318530717958679 0.217706841984230459 0.0100448605046154719
<delaware:/usr2/makino/papers/Interface>euler1
Enter k, nsteps :1 1000
k, period, nsteps, dt = 1 6.28319 1000 0.00628319
0.00628318530717958661 0.0199349143076453517 8.4329693743111499e-05
<delaware:/usr2/makino/papers/Interface>euler1
Enter k, nsteps :1 10000
k, period, nsteps, dt = 1 6.28319 10000 0.000628319
0.000628318530717958618 0.00197586995377499117 8.28467565412249203e-07
<delaware:/usr2/makino/papers/Interface>euler1
Enter k, nsteps :1 100000
k, period, nsteps, dt = 1 6.28319 100000 6.28319e-05
6.28318530717958564e-05 0.000197411570732156289 8.26997353232364789e-09
<delaware:/usr2/makino/papers/Interface>

$x$ のほうはステップを 10 倍にすると大体 1/10 になってますが、 $v$ は2桁づつ減ってます。これで合っているんでしょうか？

赤木： 「合ってるんでしょうか」って、あなたはどう思うの？

学生B：え、そう言われても、、、オイラー法は一次精度ですよね。一次精度 ということは、同じ時間積分した時の誤差がステップサイズに比例するという ことだから、少なくとも $x$ のほうはそうなってるし、それでいいと思います。

赤木： まあ、そうなんだけど、折角線型の問題を持って来たから、も うちょっとやってみて。$(x,v)$ をベクトル $\mathbf{x}$ と書くと、 $\mathbf{x}_i$ から $\mathbf{x}_{i+1}$ への変換は線型変換だから行列で書けるわね。面倒だから $k=1$ にしましょ、周期で時間と速度をスケールすれば同じことだから。

学生B：ええと、

$$\mathbf{x}_{i+1} = \pmatrix{ 1 & \Delta t \cr -\Delta t & 1 } \mathbf{x}_i$$ (6)

ですね。

赤木： そうそう、で、これ回転と定数倍に分かれるでしょ？

学生B：あ、そうですね

$$\pmatrix{ 1 & \Delta t \cr -\Delta t & 1 } = \sqrt{1+\D... ... \cos \theta & \sin \theta \cr -\sin \theta & \cos \theta}$$ (7)

で、$\theta$ は $\tan \theta = \Delta t$ の解です。あ、これなら数値解 がどうなるか分かりますね、$n$ステップ後の解は初期条件の $(1,0)$ に上の 行列の $n$ 乗を掛けるだけだから、 $(1+\Delta t^2)^{n/2}(\cos n \theta, \sin n \theta)$ ですね。

赤木： そうね、後、 $n\Delta t = 2\pi$ として、それから $\tan^{-1}$ のテイラー展開

$$\tan^{-1} x = x - x^3/3 + x^5/5 \cdots$$ (8)

を使って展開して、 $(1+\Delta t^2)^{n/2}$も展開して最初の項だけ残してみて。

学生B：ええと、、、、 $(1+2\pi^2/n, 8\pi^3/3n^2)$ となりますね。これと 上の数値解を比べると、ああ、 $n$ が大きいとよく合いますね。

赤木： そういうことね。調和振動をオイラー法で解く時には、真の解は 回転なわけで数値解の誤差は回転角の誤差とベクトルの長さというか絶対値の誤差になるわけ ね。1ステップでの回転角の誤差は $\Delta t^3/3$ で、絶対値の誤差は $\Delta t^2/2$ だから、1周期後は絶対値は $1/n$、角度は $1/n^2$ に比例 した誤差が見えるわけ。で、これがそれぞれ $x$ と $v$ の誤差になってるの よ。

学生B：そうですね、でも、結局全部手計算で答がわかるんならプログラム書かなくてもよかったんでは、、、

赤木： まあ、この問題についてはもちろんそうなわけ。数値積分するの は、解析的に解けない問題をやるためだから。でも、実際あなたの最初のプロ グラムは解析的に出せる答も出せてなかったわけよね。で、解析解と誤差の大 きさまで合ってればプログラムが正しいのはさすがに確かでしょ。

学生B：それはまあそうなんですが、、、