東京理科大学 infoserv［更新日］2000.5.1

問１３　次の Fortran プログラムの説明及びプログラムを 読んで，設問 1，2 に答えよ。

〔プログラムの説明〕

区間［a，b］における関数 f (x) の定積分 T を，ロンバーグ積分法に よって数値的に求めるプログラムである。
(1) 　 整数値 n が与えられたとき，区間［a，b］を k (＝2n) 等分して分点を x0 (＝a)，x1，
x2，…，xk−1，xk(＝b) とする。 各分点に対する関数値 f (x) を基に， 積分値を台形の面積によって近似する。例えば，区間［x0，x1］に対しては， 
                     h                                 b-a
        s(x0，x1) = ---- ｛f (x0) ＋ f (x1)｝，    h = ----
                     2                                  k

と近似する。これを区間［a，b］にわたり，すべて足し合わせた

Tn,0 ＝ s(x0，x1) ＋ s(x1，x2) ＋ … ＋ s(xk−2，xk−1) ＋ s(xk−1，xk)
が，区間を 2n 等分したときの積分値の第 0 近似である。 n ＝ 0，1，2，… として繰り返し計算して収束解を求める。

(2) 　 更に次の漸化式（リチャードソン補外）を用いることで，計算をより速く収束させる。

Tn,m ＝ (4m Tn＋1, m−1 − Tn, m−1) ／ (4m − 1), (1≦m≦n)
　例えば，n ＝ 0，1，2 の 3 通りについて，台形近似の積分値 T0,0，T1,0，T2,0を得た場合， 漸化式に従い T0,0 と T1,0 から T0,1 が， T1,0 と T2,0 から T1,1 が得られる。 更に，得られた T0,1 と T1,1を用いて， T0,2を得ることができる (図１参照)。

(3) 　 プログラムには，n の最大値 NMAX と近似精度 EPS とが与えられている。 収束条件として，|T0,n−1 − T0,n| が EPS よりも小さくなったときに， そのときの T0,nを積分の近似値として表示し終了する。n が NMAX を超えたときには， 警告メッセージを表示して終了する。
(4) 　 関数 f (x) は，外部副プログラムとして与えられる。


図１　ロンバーグ積分

〔プログラム〕 

（行番号）
01        program Romberg
02        integer, parameter :: NMAX=4, DOUBLE=kind(1.0d0)
03        real (kind=DOUBLE), parameter :: EPS=1.0d-10
04        real (kind=DOUBLE) :: a,b,h,f,fab,p,x,sum,del
05        real (kind=DOUBLE), dimension(0:NMAX,0:NMAX) :: t
06        integer :: i,k,n,m
07 !
08        a=0.0
09        b=1.0
10        del=1.0
11        h=b-a
12        fab=(f(a)+f(b))/2
13        t(0,0)=h*fab
14        write(*,'(i5,f20.15)') 0,t(0,0)
15 !
16        n=0
17        k=1
18        do while (del >= EPS)
19          n=n+1
20          k=k*2
21          h=h*0.5
22          if(n > NMAX) then
23            write(*,*) ' no convergence'
24            stop
25          endif
26          x=a+h
27          sum=fab
28          do i=1,k-1
29            sum=sum+f(x)
30            x=x+h
31          enddo
32          t(n,0)=sum*h
33          do m=1,n
34            p=4.0**m
35            t(n-m,m)=(p*t(n-m+1,m-1)-t(n-m,m-1))/(p-1.0)
36          enddo
37          write(*,'(i5,f20.15)') n,t(0,n)
38          del=abs(t(0,n)-t(0,n-1))
39        enddo
40        write(*,100) t(0,n),del
41   100  format(' integral:',f20.15,' error:',f20.15)
42        end program Romberg
43 !
44        function f(x)
45        real (kind=kind(1.0d0)) :: f
46        real (kind=kind(1.0d0)), intent(in) :: x
47        f=1.0d0/(1.0d0+x*x)
48        end function f
設問１ 　 プログラム 35 行目において，n＝2，m＝1 のときに 計算される近似式として正しい答えを，解答群の中から選べ。

解答群

ア　T0,1 ＝ (4 T1,0 − T0,0) ／ 3

イ　T1,1 ＝ (4 T2,0 − T1,0) ／ 3

ウ　T0,1 ＝ (16 T1,0 − T0,0) ／ 15

エ　T1,1 ＝ (16 T2,0 − T1,0) ／ 15

設問２ 　 次の記述中の に入れる正しい答えを，解答群の中から選べ。

このプログラムを実行したところ，図 2 のように“no convergence”と いうメッセージを出力して終了した。このとき，関数副プログラム f (x) は 回呼び出された。

プログラムでは台形近似の分点数が増えたとき， すべての分点について関数値を計算している。 しかし，そのうちの半分は既に計算済みなので，それを利用するようにしたい。 そのためには 27〜32 行目をプログラム に変更する。これによって出力結果は変わらず，関数 f (x) は 回呼び出される。

　　0　0.750000000000000　　
　　1　0.783333333333333　　
　　2　0.785529411764706　　
　　3　0.785396445940468　　
　　4　0.785398166319429　　
　　no convergence

図２　プログラムの出力例

a，c に対する解答群

ア　11　　　　イ　13　　　　ウ　15　　　　エ　17

オ　26　　　　カ　28　　　　キ　30　　　　ク　32 

b に対する解答群

ア　 27   sum=0.0
     28   do i=1,k-1,2
     29     sum=sum+f(x)
     30     x=x+h
     31   enddo
     32   t(n,0)=t(n-1,0)/2+sum*h

イ　 27   sum=fab
     28   do i=1,k-1,2
     29     sum=sum+f(x)
     30     x=x+h
     31   enddo
     32   t(n,0)=t(n-1,0)/2+sum*h

ウ　 27   sum=0.0
     28   do i=1,k-1,2
     29     sum=sum+f(x)
     30     x=x+2*h
     31   enddo
     32   t(n,0)=t(n-1,0)/2+sum*h

エ　 27   sum=fab
     28   do i=1,k-1,2
     29     sum=sum+f(x)
     30     x=x+2*h
     31   enddo
     32   t(n,0)=t(n-1,0)/2+sum*h
東京理科大学 infoserv 戻る 次頁：問14