情報地球科学演習 (2024年前学期) 第12回

今回のテーマ: 常微分方程式の初期値問題の数値解法

演習の後半戦では、物理現象の簡単な数値シミュレーションの入口にあたる実習を行っている。 その主旨は、 手計算でやるのが大変な問題を、コンピュータに代わりにやってもらう ことであるともいえる。 これまで、非線形方程式の数値解法、及び関数の定積分の計算を扱ってきたが、今回は 常微分方程式の初期値問題を、コンピュータに代わりに解いてもらうことを考える。

例題として、

• \(x\) に関する常微分方程式 \(\dfrac{d^2{x}}{d{t^2}}=-1\) を、
• \(0\le{t}\le2\) の範囲で
• \(t=0\) のとき \(x=0\) かつ \(\dfrac{d{x}}{d{t}}=1\) という初期条件のもとで

解くことを考える。 なおこの問題は、重力のある中で投げ上げられた物体の軌跡 を求める問題である。

これは2階の常微分方程式であり、講義で紹介した方法をそのまま使うことはできない。 しかし \(v\equiv\dfrac{d{x}}{d{t}}\) という新たな変数を導入してやると、

• \(v\) に関する1階の常微分方程式 \(\dfrac{d{v}}{d{t}}=-1\) かつ \(t=0\) のとき \(v=1\)
• \(x\) に関する1階の常微分方程式 \(\dfrac{d{x}}{d{t}}=v\) かつ \(t=0\) のとき \(x=0\)

という2つの問題に分割することができる。 なお、ここで導入した \(v\) はもちろん物体の速度を表している。 また、これらの方程式の解はもちろん \(v=1-t\)、\(x=t-\dfrac{1}{2}t^2\) となる。

オイラー法による解法 (その1)

まず最初は、速度 v に関する1階の常微分方程式 (\(\dfrac{dv}{dt}=-1\) かつ \(v(t=0)=1\)) だけを、オイラー法により解くことを考える。

必要な手順 (の要点) を以下に示す。 第11回 に扱った、定積分を計算する手順とよく似ていることに注意。

• オイラー法の計算に必要な「器」を宣言する。 例えば「小区間」の幅を表わす dt、及び積分区間の最小値を表わす t_kaishi と最大値を表わす t_shuryo など。 さらにここでは、t_kaishi と t_shuryo は 定数 (パラメータ) として宣言 し、同時にその値を設定 している。 parameter という「属性」の詳しい説明は ここ を参照のこと。

• 「小区間」の数や、端点での t や v の値を入れる「器」を宣言する。 ここでは、

  • 「小区間」の数 を表わす整数型変数 kosuu を 定数 (パラメータ) として宣言 し、同時にその値を設定 している。

  • 「小区間」の端点での t や v の値 を入れるための 実数型配列 を宣言し、同時に定数 kosuu を使って配列の大きさを指定 している。 なお、ここでも配列の添字の範囲の指定を (0:kosuu) とすることにより、その配列の要素は 0番めから kosuu 番め として参照され、要素の数は kosuu+1 個となる。 配列の宣言のしかたの詳しい説明は ここ を参照のこと。

  としている。

• オイラー法の計算に必要な変数に値を入れる。 ここでは「小区間」の幅 dt、端点での t の値、及び v の初期値に値を設定している。 当然ながら、「小区間」の数が必要な時は定数 kosuu を使って指定 するようにすること。

• オイラー法により、v(1) からv(kosuu) までに正しい値を計算して入れる。 講義のノートに基づいて、正しい手順を考えてプログラムに書き入れよ。 この際、dt だけ時間が進む間に、v(i) から v(i+1) へ変化する ということ、あるいは (ほとんど答なのだが) \[-1=\left[\dfrac{dv}{dt}\right]_{(i)}\simeq\frac{v_{(i+1)}-v_{(i)}}{dt}\] という関係式を意識すると理解しやすい (と思う)。 またここでも当然、配列の要素の数が必要な時は定数 kosuu を使って指定 するようにすること。

• 計算結果を表示する。 ここでは v の値だけでなく、対応する t の値も合わせて書き出すことにする。

p

r

o

g

r

a

m

e

n

s

h

u

1

3

a

i

m

p

l

i

c

i

t

n

o

n

e

i

n

t

e

g

e

r

,

p

a

r

a

m

e

t

e

r

:

:

k

o

s

u

u

=

1

6

r

e

a

l

,

d

i

m

e

n

s

i

o

n

(

0

:

k

o

s

u

u

)

:

:

t

,

v

r

e

a

l

:

:

d

t

r

e

a

l

,

p

a

r

a

m

e

t

e

r

:

:

t

_

k

a

i

s

h

i

=

0

.

0

,

t

_

s

h

u

r

y

o

=

2

.

0

i

n

t

e

g

e

r

:

:

i

d

t

=

(

t

_

s

h

u

r

y

o

-

t

_

k

a

i

s

h

i

)

/

r

e

a

l

(

k

o

s

u

u

)

d

o

i

=

0

,

k

o

s

u

u

t

(

i

)

=

t

_

k

a

i

s

h

i

+

d

t

*

r

e

a

l

(

i

)

e

n

d

d

o

v

(

0

)

=

1

.

0

d

o

i

=

0

,

k

o

s

u

u

-

1

v

(

i

+

1

)

=

e

n

d

d

o

d

o

i

=

0

,

k

o

s

u

u

w

r

i

t

e

(

*

,

*

)

t

(

i

)

,

v

(

i

)

e

n

d

d

o

e

n

d

p

r

o

g

r

a

m

e

n

s

h

u

1

3

a

これを基にして、オイラー法により上記の常微分方程式を解くプログラムを作成せよ [みほん]。 このプログラムを (例えば) enshu13a.f90 という名前で保存し、例の如く gfortran enshu13a.f90 とした後に、

$ ./a.out

と実行せよ。 また期待通りの答が得られているかを、以下の手順により確認せよ。

• 出力データの先頭3行はどうなっているか?

  $ ./a.out | head -n 3

  としたときの出力は以下のようになっているはずである。

     0.00000000       1.00000000    
    0.125000000      0.875000000    
    0.250000000      0.750000000    
        

• 出力データの末尾3行はどうなっているか?

  $ ./a.out | tail -n 3

  としたときの出力は以下のようになっているはずである。

     1.75000000     -0.750000000    
     1.87500000     -0.875000000    
     2.00000000      -1.00000000    
        

オイラー法による解法 (その2)

次は、上で作ったプログラムを改良し、速度 v に加えて位置 x も解く機能を追加 してみよう。 それには \(x\) に関する常微分方程式 \(\dfrac{d{x}}{d{t}}=v(t)\) かつ \(x(t=0)=0\) も解くようにすればよいのだが、そのために新たに追加すべき手順 (の要点) を以下に示す。

• 「小区間」の端点での x の値 を入れるための 実数型配列 を宣言する。

  r

  e

  a

  l

  ,

  d

  i

  m

  e

  n

  s

  i

  o

  n

  (

  0

  :

  k

  o

  s

  u

  u

  )

  :

  :

  x

• オイラー法の計算に必要な変数に値を入れる。 ここでは「小区間」の端点での x の初期値に値を設定している。

  		v	(	0	)		=		1	.	0
  		x	(	0	)		=		0	.	0

• オイラー法により、x(1) から x(kosuu) までに正しい値を計算して入れる。 講義のノートに基づいて、正しい手順を考えてプログラムに書き入れよ。 またここでも、dt だけ時間が進む間に、x(i) から x(i+1) へ変化する ということ、あるいは (ほとんど答なのだが) \[v_{(i)}=\left[\dfrac{dx}{dt}\right]_{(i)}\simeq\frac{x_{(i+1)}-x_{(i)}}{dt}\] という関係式を意識すると理解しやすい (と思う)。

  d

  o

  i

  =

  0

  ,

  k

  o

  s

  u

  u

  -

  1

  v

  (

  i

  +

  1

  )

  =

  x

  (

  i

  +

  1

  )

  =

  e

  n

  d

  d

  o

• 計算結果を表示する。

  d

  o

  i

  =

  0

  ,

  k

  o

  s

  u

  u

  w

  r

  i

  t

  e

  (

  *

  ,

  *

  )

  t

  (

  i

  )

  ,

  v

  (

  i

  )

  ,

  x

  (

  i

  )

  e

  n

  d

  d

  o

これらを参考にして、先に作ったプログラム enshu13a.f90 に、x を解く機能を追加したプログラムを作成せよ [みほん]。 そのプログラムを用いて (例の如く gfortran enshu13b.f90 などとして a.out を新たに作った後に)

$ ./a.out

と実行し、期待通りの答が得られているかどうかを確認せよ。 プログラムが正しく作られていれば、出力の先頭3行は以下のようになっているはずである。

   0.00000000       1.00000000       0.00000000    
  0.125000000      0.875000000      0.125000000    
  0.250000000      0.750000000      0.234375000    

正しく動いていることが確認できたら、本日の演習の出席確認のため、このプログラムを 末尾のフォーム 経由で提出せよ。

正しい出力が得られていることが確認できたら、

$ ./a.out > tvx.txt

と実行することにより、出力を tvx.txt という名前のファイルに書き出させておくこと。

オイラー法による解法 (その3)

上で求めた数値的な答がどれだけ正確なものかを、Gnuplotでグラフを描くことにより調べてみよう。

Gnuplot を起動し、まずは以下のように入力してみよ。

gnuplot> plot "tvx.txt" using 1:2

とすると、横軸に t、縦軸に v をとったグラフが表示される。 さらに

gnuplot> plot "tvx.txt" using 1:2,1-x

とすることで、v の真の答と合わせたグラフが表示される。 両者の比較により、数値的に求めた v がどの程度正確であるかを確かめよ。

同様に

gnuplot> plot "tvx.txt" using 1:3

とすると、横軸に t、縦軸に x をとったグラフが表示される。 また

gnuplot> plot "tvx.txt" using 1:3,x-x**2/2

とすることで、x の真の答と合わせたグラフが表示される。 両者の比較により、数値的に求めた x がどの程度正確であるかを確かめよ。

オイラー法による解法 (その4)

上で作ったプログラムで得られた x の計算結果は、真の答にさほど近いものではなかったはずだ。 実際、t=2 における真の答は x=0 であるのに対し、数値的な答では x は 0 から大きく外れていた。 当然予想されるように、「小区間」の幅を小さくすれば、言い変えれば 「小区間」の数を多くすれば、真の答により近い答が出てくる はずである。 そこで、「小区間」の数と、計算値の正確さの間の関係 を確かめてみよう。

「小区間」の数 kosuu を 2 から 4096 あるいは 8192 程度まで2倍ずつ増やしながら同じ計算 を行い、数値的に求めた x(t=2) と真の答からのずれがどう変化しているか を調べてみよ [みほん]。 なおそのためには、

$ ./a.out | tail -n 1

のようにプログラムを実行し、出力の最後の1行のみを表示させてやると便利である。 さらに kosuu と真の答からのずれの間にはどのような関係があるか。 必要に応じて両対数グラフにプロットしてみるなどして、その傾向を確認せよ。

ここで見たように、オイラー法では、「小区間」の数を増やしても、なかなか計算の精度が上がってくれない。 このことは、オイラー法で精度よく計算するためには、「小区間」を非常に細かくとる必要がある ことを意味している。 という訳でオイラー法は (「常微分方程式の解法の入門」という意味はあるものの) 実用的とはいえない。 実際の問題では、オイラー法よりも数段「高精度な」方法、例えば ルンゲ-クッタ法 と呼ばれる一連の方法を用いることがほとんどである。

余力のある人向けのおまけ: オイラー法による解法 (その5)

次回の演習では、オイラー法とその他の高精度な方法との間で、得られる答の正確さがどれくらい違うかを検討していく予定にしている。 そのための (ちょっと複雑な) プログラムを作る上では、微分方程式の右辺の式を、自前の関数として定義しておく のが便利である。 そこで (来週の課題を一部先取りするつもりで) 余力のある人は以下のプログラム例を参考にして、(その2) で作ったプログラムを修正し、自前の関数を定義・使用する仕様で書き直してみよう。 ただし (例によって) 肝心なところは (わざと) 抜いているので、各自で考えて埋めてほしい。

なお、このプログラムのミソは以下の通り。

• contains 以下で、v を時間 dt だけ進めるのに使う自前の関数 dvdt を宣言し、しかるべく定義する。 この関数 dvdt で、\(\dfrac{d{v}}{d{t}}=\)何々 の式の右辺に相当する式を計算させる。 ただし今回の例では、右辺の計算に x が出てこないのだが、次回の課題への流用を意図して x を関数の引数に含めてある。

  同様に、x を時間 dt だけ進めるのに使う自前の関数 dxdt を宣言し、\(\dfrac{d{x}}{d{t}}=\)何々 の式の右辺に相当する式を計算するよう定義する。

• しかるべく定義された dvdt や dxdt を使って、v(1) からv(kosuu)、x(1) から x(kosuu) の値を計算させる。

また、このような書き方をすることで、オイラー法の前提である以下の近似的な関係式 \[ \dfrac{d{v}}{d{t}}\simeq\dfrac{v_{(i+1)}-v_{(i)}}{dt} ,\quad \dfrac{d{x}}{d{t}}\simeq\dfrac{x_{(i+1)}-x_{(i)}}{dt} \] がより明確に意識できるようになるだろう (と思う)。

p

r

o

g

r

a

m

e

n

s

h

u

1

3

b

f

i

m

p

l

i

c

i

t

n

o

n

e

i

n

t

e

g

e

r

,

p

a

r

a

m

e

t

e

r

:

:

k

o

s

u

u

=

1

6

r

e

a

l

,

d

i

m

e

n

s

i

o

n

(

0

:

k

o

s

u

u

)

:

:

t

,

v

,

x

r

e

a

l

:

:

d

t

r

e

a

l

,

p

a

r

a

m

e

t

e

r

:

:

t

_

k

a

i

s

h

i

=

0

.

0

,

t

_

s

h

u

r

y

o

=

2

.

0

i

n

t

e

g

e

r

:

:

i

d

t

=

(

t

_

s

h

u

r

y

o

-

t

_

k

a

i

s

h

i

)

/

r

e

a

l

(

k

o

s

u

u

)

d

o

i

=

0

,

k

o

s

u

u

t

(

i

)

=

t

_

k

a

i

s

h

i

+

d

t

*

r

e

a

l

(

i

)

e

n

d

d

o

v

(

0

)

=

1

.

0

x

(

0

)

=

0

.

0

d

o

i

=

0

,

k

o

s

u

u

-

1

v

(

i

+

1

)

=

v

(

i

)

+

d

t

*

d

v

d

t

(

x

(

i

)

)

x

(

i

+

1

)

=

x

(

i

)

+

d

t

*

d

x

d

t

(

v

(

i

)

)

e

n

d

d

o

d

o

i

=

0

,

k

o

s

u

u

w

r

i

t

e

(

*

,

*

)

t

(

i

)

,

v

(

i

)

,

x

(

i

)

e

n

d

d

o

c

o

n

t

a

i

n

s

f

u

n

c

t

i

o

n

d

v

d

t

(

x

)

r

e

a

l

,

i

n

t

e

n

t

(

i

n

)

:

:

x

r

e

a

l

:

:

d

v

d

t

d

v

d

t

=

e

n

d

f

u

n

c

t

i

o

n

d

v

d

t

f

u

n

c

t

i

o

n

d

x

d

t

(

v

)

r

e

a

l

,

i

n

t

e

n

t

(

i

n

)

:

:

v

r

e

a

l

:

:

d

x

d

t

d

x

d

t

=

v

e

n

d

f

u

n

c

t

i

o

n

d

x

d

t

e

n

d

p

r

o

g

r

a

m

e

n

s

h

u

1

3

b

f

出席確認用プログラム提出フォーム

Moodle のコースページ 以下の 課題ファイル提出 により、各自で作成した (enshu13b.f90 に相当する) オイラー法によって本日の課題の常微分方程式を解く Fortran 90 プログラムの最終版を提出せよ。