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

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

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

今日は新たな例題として、

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

解くことを考える。なおこの常微分方程式は、単振動 を記述するものである。今日の具体的な目標は、この方程式の答をオイラー法、ホイン法や中点法、さらにはルンゲ-クッタ法を使って求め、それらの正確さを比較してみることである。

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

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

という2つの問題に分割することができる。また、これらの方程式の解はもちろん $x=\sin(\pi{t})$、$v=\pi\cos(\pi{t})$ となる。

オイラー法による解法

まず最初は、オイラー法を使ってこの方程式を解いてみよう。必要な手順 (の要点) を以下に示す。基本的には第12回に作成した、オイラー法により物体の運動の軌跡を求めるプログラムの一部を修正 すればよい。なお以下では、「自前の関数」を使わずに書いたプログラムをもとにして説明しているが、「自前の関数」を使ったプログラムを使いたい人は、これにならって適切に修正してほしい。

オイラー法の計算に必要な「器」を宣言する。例えば「小区間」の幅を表わす dt、及び積分区間の最小値を表わす t_kaishi と最大値を表わす t_shuryo など。さらにここでは、t_kaishi と t_shuryo は 定数 (パラメータ) として宣言 し、同時にその値を設定 している。 parameter という「属性」の詳しい説明はここを参照のこと。またついでに、円周率 π を表わす pi、及びその2乗を表わす pi2jo という定数も宣言し、値を設定している。
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
r e a l , p a r a m e t e r : : p i = 4 . 0 * a t a n ( 1 . 0 ) , p i 2 j o = p i * * 2
「小区間」の数や、端点での t や v、x の値を入れる「器」を宣言する。ここでは、
- 「小区間」の数 を表わす整数型変数 kosuu を 定数 (パラメータ) として宣言 し、同時にその値を設定 している。
- 「小区間」の端点での t や v、x の値 を入れるための 実数型配列 を宣言し、同時に定数 kosuu を使って配列の大きさを指定 している。なお、ここでも配列の添字の範囲の指定を (0:kosuu) とすることにより、その配列の要素は 0番めから kosuu 番め として参照され、要素の数は kosuu+1 個となる。配列の宣言のしかたの詳しい説明はここを参照のこと。
としている。
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
オイラー法の計算に必要な変数に値を入れる。ここでは「小区間」の幅 dt、端点での t の値、及び v と x の初期値に値を設定している。当然ながら、「小区間」の数が必要な時は定数 kosuu を使って指定 するようにすること。
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 ) = p i
x ( 0 ) = 0 . 0
オイラー法により、v(1) から v(kosuu) まで、及び x(1) から x(kosuu) までに正しい値を計算して入れる。講義のノートに基づいて、正しい手順を考えてプログラムに書き入れよ。ここでも当然、配列の要素の数が必要な時は定数 kosuu を使って指定 するようにすること。

d o i = 0 , k o s u u - 1
v ( i + 1 ) = v ( i ) + d t *
x ( i + 1 ) = x ( i ) + d t *
e n d d o

ただし、「自前の関数」を使ったプログラムをもとにして修正する場合には、ここを修正するのではなく、contains 以下にある「自前の関数」の内容を修正するだけ でよい (これが「自前の関数」使うメリットの1つ)。
計算結果を表示する。ここでは v や x の値だけでなく、対応する t の値も合わせて書き出すことにする。
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

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

$ ./a.out

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

   0.00000000       3.14159274       0.00000000    
  0.125000000       3.14159274      0.392699093    
  0.250000000       2.65711975      0.785398185

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

$ ./a.out > tvx2.txt

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

さらに、上で求めた数値的な答がどれだけ正確なものかを、Gnuplotでグラフを描くことにより調べてみよう。 Gnuplot を起動し、まずは以下のように入力してみよ。

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

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

gnuplot> plot "tvx2.txt" using 1:2,pi*cos(pi*x)

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

次に

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

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

gnuplot> plot "tvx2.txt" using 1:3,sin(pi*x)

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

ホイン法あるいは中点法による解法

講義でも述べた通り、オイラー法は実用的な方法ではなく、実際の問題を解くにはより正確さに優れた別の方法を使うべきである。そこで、オイラー法による計算結果と、より正確な方法で得られた計算結果を比較し、答の正確さの違いを調べてみよう。

ここでは、上で作ったオイラー法によるプログラムを改良し、「2次の解法」である ホイン法あるいは中点法によって解く機能を追加 してみよう。新たに追加する手順 (の要点) を以下に示す。

ホイン法あるいは中点法で求めた v や x の値を入れる「器」を宣言する。ここではそれぞれ v2 と x2 という名前の配列を用意している。
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 i m e n s i o n ( 0 : k o s u u ) : : v 2 , x 2
ホイン法あるいは中点法の計算に必要な「器」を宣言する。以下では、これらの方法で中間的に用いる k₁ と k₂ にあたるものを宣言している。ただしこの問題では、v に関する微分方程式と x に関する微分方程式の2つを連立させて解く ことになるので、 k₁ や k₂ も v 用と x 用の2つを用意 しておく必要がある。
r e a l : : k 1 v , k 2 v , k 1 x , k 2 x

ホイン法あるいは中点法により、v2(0) から v2(kosuu) まで、及び x2(0) から x2(kosuu) までに正しい値を計算して入れる。講義のノートに基づいて、正しい手順を考えてプログラムに書き入れよ。ここでは当然、速度の値を参照したい時は配列 v2 を、同様に 位置の値を参照したい時は配列 x2 を使って指定するようにすること。また後の作業の都合上、オイラー法で答を求める手順も消さずに置いておく ことにしよう。

ホイン法を使う場合

中点法を使う場合

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

計算結果を表示する。ここでは、オイラー法で求めた結果と合わせ、ホイン法あるいは中点法で求めた結果も表示させる。
d o i = 0 , k o s u u
w r i t e ( * , * ) t ( i ) , v ( i ) , x ( i ) , v 2 ( i ) , x 2 ( i )
e n d d o

これを基にして、オイラー法に加えてホイン法あるいは中点法を用いて上記の常微分方程式を解くプログラムを作成せよ [ホイン法のみほん (関数不使用版と使用版)] [中点法のみほん (関数不使用版と使用版)]。このプログラムを (例えば) enshu14b.f90 という名前で保存し、例の如く gfortran enshu14b.f90 とした後に、

$ ./a.out

と実行し、期待通りの答が得られているかどうかを確認せよ。プログラムが正しく作られていれば、出力の先頭3行は以下のようになっているはずである。 (このうち4列目と5列目がホイン法あるいは中点法による答)。

ホイン法を使う場合

   0.00000000       3.14159274       0.00000000       3.14159274       0.00000000    
  0.125000000       3.14159274      0.392699093       2.89935613      0.392699093    
  0.250000000       2.65711975      0.785398185       2.19132447      0.724839091

中点法を使う場合

   0.00000000       3.14159274       0.00000000       3.14159274       0.00000000    
  0.125000000       3.14159274      0.392699093       2.89935613      0.392699093    
  0.250000000       2.65711975      0.785398185       2.19132447      0.724839032

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

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

$ ./a.out > tvx3.txt

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

さらに、ホイン法あるいは中点法で求めた数値的な答がオイラー法の答と比べてどれだけ正確なものかを、Gnuplot でグラフを描くことにより調べてみよう。 Gnuplot を起動し、

gnuplot> plot "tvx3.txt" using 1:2, "tvx3.txt" using 1:4, pi*cos(pi*x)

とすると、横軸に t、縦軸に v をとったグラフ上に、オイラー法による結果が赤、ホイン法あるいは中点法による結果が緑、真の答が青で表示される。同様に

gnuplot> plot "tvx3.txt" using 1:3, "tvx3.txt" using 1:5, sin(pi*x)

とすると、横軸に t、縦軸に x をとったグラフ上に、オイラー法による結果が赤、ホイン法あるいは中点法による結果が緑、真の答が青で表示される。これらのグラフを比較することにより、「小区間」への分割数が同じであっても、ホイン法や中点法ではオイラー法と比べて格段に正確さが増していることが分かるだろう。

ルンゲ-クッタ法による解法

最後に、「4次の解法」である ルンゲ-クッタ法によって解く機能も追加 したプログラムを作成してみよう。そしてオイラー法、ホイン法や中点法と結果を比較し、答がさらに正解に近くなるようすを確認してみよう。

新たに追加する手順 (の要点) を以下に示す。

ルンゲ-クッタ法で求めた v や x の値を入れる「器」を宣言する。ここではそれぞれ v4 と x4 という名前の配列を用意している。
r e a l , d i m e n s i o n ( 0 : k o s u u ) : : t , v , x , v 2 , x 2
r e a l , d i m e n s i o n ( 0 : k o s u u ) : : v 4 , x 4
ルンゲ-クッタ法の計算に必要な「器」を宣言する。以下では、ルンゲ-クッタ法で中間的に用いる k₁ から k₄ にあたるものを宣言している。ただしホイン法の場合と同様に、 k₁ から k₄ も v 用と x 用の2つを用意 しておく必要がある。
r e a l : : k 1 v , k 2 v , k 3 v , k 4 v , k 1 x , k 2 x , k 3 x , k 4 x

ルンゲ-クッタ法により、v4(0) から v4(kosuu) まで、及び x4(0) から x4(kosuu) までに正しい値を計算して入れる。講義のノートに基づいて、正しい手順を考えてプログラムに書き入れよ。ここでは当然、速度の値を参照したい時は配列 v4 を、同様に 位置の値を参照したい時は配列 x4 を使って指定するようにすること。また後の作業の都合上、オイラー法やホイン法あるいは中点法で答を求める手順も消さずに置いておく ことにしよう。

v	(	0	)		=		p	i
x	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			v	(	i	+	1	)		=
			x	(	i	+	1	)		=
e	n	d		d	o
v	2	(	0	)		=		p	i
x	2	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			k	1	v		=
			k	1	x		=
			k	2	v		=
			k	2	x		=
			v	2	(	i	+	1	)		=		v	2	(	i	)	+
			x	2	(	i	+	1	)		=		x	2	(	i	)	+
e	n	d		d	o
v	4	(	0	)		=		p	i
x	4	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			k	1	v		=
			k	1	x		=
			k	2	v		=
			k	2	x		=
			k	3	v		=
			k	3	x		=
			k	4	v		=
			k	4	x		=
			v	4	(	i	+	1	)		=		v	4	(	i	)	+	(	k	1	v	+	2	.	0	*	k	2	v	+	2	.	0	*	k	3	v	+	k	4	v	)	/	6	.	0
			x	4	(	i	+	1	)		=		x	4	(	i	)	+	(	k	1	x	+	2	.	0	*	k	2	x	+	2	.	0	*	k	3	x	+	k	4	x	)	/	6	.	0
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 ) , v 2 ( i ) , x 2 ( i ) , v 4 ( i ) , x 4 ( i )
e n d d o

上で作成したプログラムに、ルンゲ-クッタ法で同じ方程式を解く機能を追加したプログラムを作成せよ [みほん (関数不使用版と使用版)]。このプログラムを (例えば) enshu14c.f90 という名前で保存し、例の如く gfortran enshu14c.f90 とした後に、

$ ./a.out

と実行し、期待通りの答が得られているかどうかを確認せよ。プログラムが正しく作られていれば、出力の先頭3行は以下のようになっているはずである (このうち6列目と7列目がルンゲ-クッタ法による答)。

   0.00000000       3.14159274       0.00000000       3.14159274       0.00000000       3.14159274       0.00000000    
  0.125000000       3.14159274      0.392699093       2.89935613      0.392699093       2.90246916      0.382605910    
  0.250000000       2.65711975      0.785398185       2.19132447      0.724839091       2.22165728      0.706967354

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

$ ./a.out > tvx4.txt

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

さらに、ルンゲ-クッタ法で求めた数値的な答がオイラー法やホイン法で求めたのもの比べてどれだけ正確なものかを、Gnuplot でグラフを描くことにより調べてみよう。 Gnuplot を起動し、

gnuplot> plot "tvx4.txt" using 1:2, "tvx4.txt" using 1:4, "tvx4.txt" using 1:6, pi*cos(pi*x)

とすると、横軸に t、縦軸に v をとったグラフ上に、オイラー法による結果が赤、ホイン法による結果が緑、ルンゲ-クッタ法による結果が青、真の答が紫で表示される。同様に

gnuplot> plot "tvx4.txt" using 1:3, "tvx4.txt" using 1:5,"tvx4.txt" using 1:7, sin(pi*x)

とすると、横軸に t、縦軸に x をとったグラフ上に、オイラー法による結果が赤、ホイン法による結果が緑、ルンゲ-クッタ法による結果が青、真の答が紫で表示される。オイラー法をホイン法に変えるだけで格段の違いがあったように、ルンゲ-クッタ法に変えることによってもさらにまた劇的な違いがもたらされることが分かるだろう。

余力のある人向けのおまけ (その1): ルンゲ-クッタ-ギル法

「2次の解法」にホイン法と中点法の2つのバリエーションがあった (本当は無限にある) のと同様に、上で挙げたルンゲ-クッタ法の他にも「4次の解法」が存在する。ここではその1つとして、ルンゲ-クッタ-ギル (Runge-Kutta-Gill) 法と呼ばれる方法を紹介する。

ここに挙げたプログラム (関数不使用版と使用版) は、今日の課題である単振動の問題をルンゲ-クッタ-ギル法で解くものである。標準的なルンゲ-クッタ法と同様にルンゲ-クッタ-ギル法でも、k₁ から k₄ までの途中の段階を計算し、それらの値をうまく組み合わせて解を求めていく。ただしルンゲ-クッタ-ギル法では、途中の段階の計算の際に登場する係数が「きれいな」有理数になっていないというところに「気持ちの悪さ」を感じてしまうだろう。とはいうもののこの方法でも、標準的なルンゲ-クッタ法と同じくらいに正確な解を求めることができる (そりゃ当然だ)。

ただし、実際にルンゲ-クッタ-ギル法を採用しているプログラムであっても、世の中にあるもののほとんどはこれとは違った手順で計算を進めているはずだろう。例えばここに挙げたプログラム (関数不使用版と使用版) は、「修正版」のルンゲ-クッタ-ギル法の手順に基づいて、今日の課題である単振動の問題を解くものである。これらのプログラムの中では、ルンゲ-クッタ-ギル法の手順を以下のように書き直した上で実行している。

${k^{(1)}}=f(t_i,x_i)\times{h}$,	$x^{(1)}={x}_i+\dfrac{1}{2}{k^{(1)}}$,	${Q^{(1)}}={k^{(1)}}$
${k^{(2)}}=f(t_i+\dfrac{1}{2}{h},{x^{(1)}})\times{h}$,	${x^{(2)}}={x^{(1)}}+\dfrac{2-\sqrt{2}}{2}({k^{(2)}}-{Q^{(1)}})$,	${Q^{(2)}}=\dfrac{3\sqrt{2}-4}{2}{Q^{(1)}}+(2-\sqrt{2}){k^{(2)}}$
${k^{(3)}}=f(t_i+\dfrac{1}{2}{h},{x^{(2)}})\times{h}$,	${x^{(3)}}={x^{(2)}}+\dfrac{2+\sqrt{2}}{2}({k^{(3)}}-{Q^{(2)}})$,	${Q^{(3)}}=-\dfrac{4+3\sqrt{2}}{2}{Q^{(2)}}+(2+\sqrt{2}){k^{(3)}}$
${k^{(4)}}=f(t_i+h,{x^{(3)}})\times{h}$,	${x}_{i+1}={x^{(3)}}+\dfrac{1}{6}{k^{(4)}}-\dfrac{1}{3}{Q^{(3)}}$

この「修正版」の手順の大事な特徴の1つは、計算の途中のある段階での $k^{(i)}$、$x^{(i)}$、$Q^{(i)}$ の値を求める際に、その直前の段階における値しか必要としていない ことである。そのため、途中の段階で求まった変数の値を全てプログラム中で別々に保存しておく必要がなくなり、プログラム中で使用する作業変数の数を減らすことができる。この特徴は、コンピュータに塔載されたメモリ量がものすごく小さかった太古の昔では非常に有意義であり、それゆえ長い歴史をもつプログラムの中にはルンゲ-クッタ-ギル法を用いているものも少なからず存在する。とはいうものの、大容量のメモリが積まれている最近のコンピュータを使う場合には、この点を気にする必要はまずないだろうし、「どうしてもルンゲ-クッタ-ギル法でなければ....」という理由もほとんどない (ような気がする)。そのせいもあってか、亀山の手元にいくつかある数値計算の教科書のうちで、ルンゲ-クッタ-ギル法を紹介しているものは、「1963年初版」というものすごく古い1冊だけだったりする。

余力のある人向けのおまけ (その2): ケプラー運動

精度のよい常微分方程式の数値解法の威力が最も顕著に現れる問題の1つは、太陽のまわりを公転する惑星の軌道を求めてみる ことである。なおこの運動は「ケプラー運動」と呼ばれ、有名なケプラーの法則に従う。

第1法則 (楕円軌道の法則): 惑星は、太陽をひとつの焦点とする楕円軌道上を動く。
第2法則 (面積速度一定の法則): 惑星と太陽とを結ぶ線分が単位時間に描く面積は一定である。
第3法則 (調和の法則): 惑星の公転周期の2乗は、軌道の長半径の3乗に比例する。

以下では簡単のため、太陽は原点にあり、惑星は x-y 平面上を運動しているものとしよう。時刻 $t$ における惑星の座標を $(x(t),y(t))$ と書くと、運動方程式は $\dfrac{d^2{x}}{d{t^2}}=-\dfrac{GMx(t)}{r(t)^3}$ かつ $\dfrac{d^2{y}}{d{t^2}}=-\dfrac{GMy(t)}{r(t)^3}$、あるいは惑星の位置ベクトル $\vec{x}=(x,y)$ を使ってまとめて書くと \[\dfrac{d^2\vec{x}}{d{t^2}}=-\dfrac{GM}{r(t)^3}\vec{x}(t)\] で与えられる。ここで $r(t)\equiv|\vec{x}(t)|=\sqrt{x(t)^2+y(t)^2}$ は時刻 t における太陽からの距離を表わしており、$G$ は万有引力定数、$M$ は太陽の質量を表す。また以下の例では簡単のため、$GM=1$ とみなしている。

このプログラム (コピペ可) を参考にして、ケプラー運動を計算するプログラムを作成してみよ。その際、計算にはオイラー法、ホイン法、ルンゲ-クッタ法の3通りを試してみよ。なおこのプログラムをコンパイル・実行すると kepler.txt という名前のファイルに計算結果が書き出されるが、これを Gnuplot で表示するには例えば

gnuplot> plot "kepler.txt" using 2:3 with lines, "kepler.txt" using 4:5

とすると、理論的に求まる楕円軌道が赤の実線で、オイラー法で求まる軌道が緑の点で表示される。同様に以下の手順によりそれぞれ、ホイン法、ルンゲ-クッタ法によって求まる軌道を表示させることができる。

gnuplot> plot "kepler.txt" using 2:3 with lines, "kepler.txt" using 6:7
gnuplot> plot "kepler.txt" using 2:3 with lines, "kepler.txt" using 8:9

この計算結果から、ケプラーの第1法則が成り立っていることを実感してみよ。また 離心率の大きな楕円軌道をうまく計算で求めるには慎重な取り扱い (時間刻みを非常に細かくとる、あるいは精度の高い解法を用いる) が求められていることを確認してみよ。

余力のある人向けのおまけ (その3): 偏微分方程式の初期値・境界値問題

物理現象を記述する偏微分方程式の簡単な例として、熱伝導現象を表わす偏微分方程式をとりあげよう。ある長さをもつ細い棒の中を、熱伝導により熱が伝わる現象を考える。このとき、棒の中の温度 $T$ は、位置 $x$ と時間 $t$ という2つの独立変数によって指定され、 $\rho{C_p}\dfrac{\partial{T(x,t)}}{\partial{t}}=k\dfrac{\partial^2{T(x,t)}}{\partial{x^2}}$ という偏微分方程式 (そのものズバリで「熱伝導方程式」と呼ばれる) に従う。ここで $\rho$ は棒の密度、$C_p$ は比熱、$k$ は熱伝導率と呼ばれる。また (容易に想像できるように)、この偏微分方程式を解いて温度分布 $T(x,t)$ を求めるには、

時刻 $t=0$ における棒の内部の温度分布 (初期条件)
棒の両端での温度 (「境界条件」 と呼ばれる)

の2つを指定してやる必要がある。このような問題は「初期値・境界値問題」と呼ばれており、物理現象の数値シミュレーションでは頻繁に登場する。

このプログラム (これもコピペ可) を参考にして、1次元熱伝導方程式を計算するプログラムを作成してみよ。その際、各項を計算する手順がプログラム中でどのように記述されているかを検討してみよ。また計算を実行し、計算で得られた「棒の中の温度分布の時間変化」をグラフに表示してみることにより、熱伝導によって熱が伝わるようすを実感してみよ。

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

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

\({k^{(1)}}=f(t_i,x_i)\times{h}\),	\(x^{(1)}={x}_i+\dfrac{1}{2}{k^{(1)}}\),	\({Q^{(1)}}={k^{(1)}}\)
\({k^{(2)}}=f(t_i+\dfrac{1}{2}{h},{x^{(1)}})\times{h}\),	\({x^{(2)}}={x^{(1)}}+\dfrac{2-\sqrt{2}}{2}({k^{(2)}}-{Q^{(1)}})\),	\({Q^{(2)}}=\dfrac{3\sqrt{2}-4}{2}{Q^{(1)}}+(2-\sqrt{2}){k^{(2)}}\)
\({k^{(3)}}=f(t_i+\dfrac{1}{2}{h},{x^{(2)}})\times{h}\),	\({x^{(3)}}={x^{(2)}}+\dfrac{2+\sqrt{2}}{2}({k^{(3)}}-{Q^{(2)}})\),	\({Q^{(3)}}=-\dfrac{4+3\sqrt{2}}{2}{Q^{(2)}}+(2+\sqrt{2}){k^{(3)}}\)
\({k^{(4)}}=f(t_i+h,{x^{(3)}})\times{h}\),	\({x}_{i+1}={x^{(3)}}+\dfrac{1}{6}{k^{(4)}}-\dfrac{1}{3}{Q^{(3)}}\)

v	(	0	)		=		p	i
x	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			v	(	i	+	1	)		=
			x	(	i	+	1	)		=
e	n	d		d	o
v	2	(	0	)		=		p	i
x	2	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			k	1	v		=
			k	1	x		=
			k	2	v		=
			k	2	x		=
			v	2	(	i	+	1	)		=		v	2	(	i	)	+
			x	2	(	i	+	1	)		=		x	2	(	i	)	+
e	n	d		d	o
v	4	(	0	)		=		p	i
x	4	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			k	1	v		=
			k	1	x		=
			k	2	v		=
			k	2	x		=
			k	3	v		=
			k	3	x		=
			k	4	v		=
			k	4	x		=
			v	4	(	i	+	1	)		=		v	4	(	i	)	+	(	k	1	v	+	2	.	0	*	k	2	v	+	2	.	0	*	k	3	v	+	k	4	v	)	/	6	.	0
			x	4	(	i	+	1	)		=		x	4	(	i	)	+	(	k	1	x	+	2	.	0	*	k	2	x	+	2	.	0	*	k	3	x	+	k	4	x	)	/	6	.	0
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	)	,	v	2	(	i	)	,	x	2	(	i	)	,	v	4	(	i	)	,	x	4	(	i	)
e	n	d		d	o

v	(	0	)		=		p	i
x	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			v	(	i	+	1	)		=
			x	(	i	+	1	)		=
e	n	d		d	o
v	2	(	0	)		=		p	i
x	2	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			k	1	v		=
			k	1	x		=
			k	2	v		=
			k	2	x		=
			v	2	(	i	+	1	)		=		v	2	(	i	)	+
			x	2	(	i	+	1	)		=		x	2	(	i	)	+
e	n	d		d	o
v	4	(	0	)		=		p	i
x	4	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			k	1	v		=
			k	1	x		=
			k	2	v		=
			k	2	x		=
			k	3	v		=
			k	3	x		=
			k	4	v		=
			k	4	x		=
			v	4	(	i	+	1	)		=		v	4	(	i	)	+	(	k	1	v	+	2	.	0	*	k	2	v	+	2	.	0	*	k	3	v	+	k	4	v	)	/	6	.	0
			x	4	(	i	+	1	)		=		x	4	(	i	)	+	(	k	1	x	+	2	.	0	*	k	2	x	+	2	.	0	*	k	3	x	+	k	4	x	)	/	6	.	0
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	)	,	v	2	(	i	)	,	x	2	(	i	)	,	v	4	(	i	)	,	x	4	(	i	)
e	n	d		d	o

v	(	0	)		=		p	i
x	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			v	(	i	+	1	)		=
			x	(	i	+	1	)		=
e	n	d		d	o
v	2	(	0	)		=		p	i
x	2	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			k	1	v		=
			k	1	x		=
			k	2	v		=
			k	2	x		=
			v	2	(	i	+	1	)		=		v	2	(	i	)	+
			x	2	(	i	+	1	)		=		x	2	(	i	)	+
e	n	d		d	o
v	4	(	0	)		=		p	i
x	4	(	0	)		=		0	.	0
d	o		i		=		0	,	k	o	s	u	u	-	1
			k	1	v		=
			k	1	x		=
			k	2	v		=
			k	2	x		=
			k	3	v		=
			k	3	x		=
			k	4	v		=
			k	4	x		=
			v	4	(	i	+	1	)		=		v	4	(	i	)	+	(	k	1	v	+	2	.	0	*	k	2	v	+	2	.	0	*	k	3	v	+	k	4	v	)	/	6	.	0
			x	4	(	i	+	1	)		=		x	4	(	i	)	+	(	k	1	x	+	2	.	0	*	k	2	x	+	2	.	0	*	k	3	x	+	k	4	x	)	/	6	.	0
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	)	,	v	2	(	i	)	,	x	2	(	i	)	,	v	4	(	i	)	,	x	4	(	i	)
e	n	d		d	o