レッスン 3

問題3-1.

f(x) = ax²+bx+cをd〜eまで積分し、分割数nの違いによる答えの変化をグラフ化するプログラムを作成しましょう。積分の方法は、台形公式でもシンプソン法でも何でも構いません。但し、a〜e、nは別のファイルから読み込むこととします。ただの積分でも良いのですが、分割数による違いを付けると２重ループになるので、付けました。いきなりで難しい人は、まずは分割数なしでやってみても構いません。

これだけだとあまりに不親切っぽいので、台形公式の考え方について説明します。積分というのは、高校生の頃習ったかと思いますが、積分する範囲をdxごとの短冊状に区切って、それを足し合わせる。で、dxを無限小にするというヤツでした。

ただ、数値計算においては、この無限小ってヤツが出来ません。数値計算の情報は全て離散量であり、厳密な連続量を表現出来ません（どんなに頑張っても点のデータで線にはなれない）。

で、しょうがないから、なるべく誤差の少ない近似を考えます。台形近似はかなり簡単なものですが、長方形近似よりは誤差が少なかろうというものです。イメージは以下の通りです。シンプソン法は自分で調べてください。

[ 図-3.1 ]

fとかの関数の設定法が分からない人はこちらを参照してください。
keyword[parameter,dimension,do]

解答例

プログラム本体 - ex3.f -

01	c
02	c	sample program for ex3.f
03	c	integral
04	c
05			parameter(ni=120)
06			dimension x(ni),f(ni)
07	c
08			open(10,file='ans.d',status='unknown')
09			open(11,file='ivalue.d',status='old')
10			read(11,*) a,b,c
11			read(11,*) d,e
12			read(11,*) nj
13			close(11)
14	c
15			do j=1,nj+1,2
16			n=j
17	c
18			dx=(e-d)/float(n)
19	c
20			do i=1,n+1
21			x(i)=d+(i-1)*dx
22			f(i)=a*x(i)**2+b*x(i)+c
23			end do
24	c
25			s=0.
26			do i=1,n
27			ds=(f(i)+f(i+1))*dx*.5
28			s=s+ds
29			end do
30	c
31			write(10,*) n,s
32	c
33			end do
34			close(10)
35			end

初期設定ファイル - ivalue.d -

01	1,2,4
02	-10,3
03	100
04
05			read(11,*) a,b,c
06			read(11,*) d,e
07			read(11,*) nj

グラフ

解説

それほど解説することがないのですが、doループのiは台形を足していくときの繰返しを表し、jで分割数を変化さています。

18行目のfloatは整数nを実数に変化させています。dx(実数)=(e(実数)-d(実数))/n(整数)とやるとおかしくなるコンパイラがあるためです。逆に実数を整数にしたい場合は、intとします。ex) int(a) 詳しくはこちらを参照。

ivalue.dと言うファイルが初期値を与えるファイルです。最初の3行だけで構いません。

GraphはSma4Winで作りました。縦軸が答え、横軸が分割数になっています。意外とあっさり収束することが分かります。