FORTRANをご存じの方がいらっしゃいましたら、間違いを指摘していただきたいです。

大学で現在FORTRANのシンプソン公式の課題が出ています。
自分で考えて下記プログラムを作ったのですが、どうしても
正しい答えが出力されません。

正しい答えとしては、分割数４とか５とかで
πに限りなく近い値になるそうです。
（下に添付した物は、分割数を4にしています。）

もしFORTRANをご存じの方がいらっしゃいましたら
どこがどのように間違っているのか、ご指摘いただけないでしょうか？
よろしくお願いいたします。

問題

4/(1＋x＾2)を、積分範囲：0～1
を、シンプソン公式を使って求めよ。

区間分割数については、各自適当に与え
分割数が多くなると精度が良くなる事を確認せよ。

シンプソン公式を使うと台形公式を使う場合より
少ない分割数で解（＝π）に近づくか確認せよ。





＊＊＊ここから＊＊＊

c numerical value integral

c trapezoid formula

c *** main program ***

real a,b
parameter (n=4)

a=0
b=1

call simpson(a,b,n,sumf)

write(*,*) 'sumf=',sumf

end

c *** end of main program ***

c *** sub program ***

subroutine simpson(a,b,n,sumf)

dimension f(1000)

m=2*n

h=(b-a)/float(m)

do 10 i=1,m+1
x=a+h*float(i-1)
f(i)=func(x)
10 continue

sumf=h/3.0*(f(1)+f(m+1)+4*f(m))

do 20 i=2,m
sumf=sumf+4.0*(h/3)*f(2*i+1)
20 continue

do 30 i=2,m
sumf=sumf+2.0*(h/3)*f(2*i)
30 continue

end

c *** end of sub program ***

c *** function sub program ***

function func(x)

func=4.0/(1.0+x**2)

return
end

c *** end of function sub program ***


＊＊＊ここまで＊＊＊


もしお分かりになる方がいらっしゃいましたら、
お手数ですがご指導いただけないでしょうか？

お願いいたします。

No.3 ベストアンサー 回答者： cametan_42 回答日時： 2007/11/30 18:16

ええと、一応念のためにシンプソン公式の擬似コードも挙げておきます。

1.関数simpsonを引数を積分下限a、積分上限b、分割数n、(可能であれば)引用する関数名funcとして定義する。
2.局所変数でm＝n/2を定義する。
3.同じく、局所変数で微小区間d＝(b－a)/nを定義する。
4.yの初期値f(a)＋f(b)を計算する。
5.一回目の繰り返し計算をj＝0～j＝m-1の間で行う。なお、yの初期値は4で出したyとする。
　i.局所変数としてx ＝ a*(2*j＋1)*dを計算する。
　ii.y_k ＝ y_{k-1}＋4*func(x)を計算する。
6.二回目の繰り返し経産をi＝1～i＝m－1の間で行う。なお、yの初期値は5.で出したyとする。
　i.局所変数としてx ＝ a＋2*j*dを計算する。
　ii.y_n ＝ y_{n-1}＋2*func(x)を計算する。
7.6.で出たyにd/3を掛けたモノが計算結果となる。

以上です。
ご自分の書いたソースコードと良く見比べてみてください。
1回目の繰り返し計算と2回目の繰り返し計算とは「ループさせる」範囲が異なる事に注意して下さい(2回目の方が一回分少ないです)。

この回答へのお礼

親切なご回答をありがとうございます。
Schemeというものを全く知らないので、
下のご指導にはちょっと焦ったのですが、
このご回答なら大丈夫そうです。
参考にさせていただきまして、
もう一度考え直して見たいと思っています。
本当にありがとうございました。

お礼日時：2007/12/01 01:31

No.4 回答者： Tacosan 回答日時： 2007/11/30 18:43

#1 です.

おお, 確かに 2つのループの範囲は変えないとダメですねぇ＞#3.
ということで適宜修正してください.

No.2 回答者： cametan_42 回答日時： 2007/11/30 17:42

fortranは良く知らないので、Schemeで書いてみました。

(define (simpson func a b n)
(let ((m (/ n 2))
(d (/ (- b a) n)))
(let loop0 ((j 0) (y0 (+ (func a) (func b))))
(if (= j m)
(let loop1 ((i 1) (y1 y0))
(if (= i m)
(exact->inexact (/ (* y1 d) 3))
(loop1 (+ i 1) (let ((x (+ a (* 2 i d))))
(+ y1 (* 2 (func x)))))))
(loop0 (+ j 1) (let ((x (+ a (* d (+ 1 (* 2 j))))))
(+ y0 (* 4 (func x)))))))))

(define (func x)
(/ 4 (+ 1 (* x x))))

これを実行すると、

> (simpson func 0 1 4)
3.1415686274509804

と確かに4分割でπに近い値が出ます。

まずはシンプソンの公式ですが、

U＝d/3*{f(a)＋f(b)＋4*Σ_{j=0}^{m-1} f(x_{2j+1})＋2*Σ_{j=1}^{m-1}f(x_{2j})}

です。(LaTeX表記みたいになってゴメンなさい)
さて、exboyさんが書いたソースを見てみると、

m=2*n

となっていますが、これがそもそも

m=n/2

の間違いなんじゃないかな、と思います。
そして、最初の繰り返し計算に入る前、exboyさんの記法で行くと、sumfの初期値は単にfunc(a)+func(b)だと思います。中で3回も繰り返し計算をしているように見えるんですが、多分2回で構いません。
そして、1回目の繰り返し計算(i＝0からi＝m-1の間)が終了してから、それを2回目の繰り返し計算のsumfの初期値として渡します。
あとは上に挙げた公式に従って、それぞれの繰り返し計算のxとsumfの定義式をしっかりと書けば動くと思いますよ。
上に挙げたSchemeのソースコードも見て、一回書き換えて試してみてください。

No.1 回答者： Tacosan 回答日時： 2007/11/29 20:50

ここで使っている変数名で Simpson の公式を書いてみましたか?

さしあたり気付いたところ:
１．sumf=h/3.0*(f(1)+f(m+1)+4*f(m))
は
sumf=h/3.0*(f(1)+f(m+1))
じゃないですか?
２．sumf を計算する 2個の do はいずれも, 範囲を 2, m ではなく 1, n にしないといけないような.

この回答へのお礼

親切におしえていただきましてありがとうございます。
4*f(m)も、最初はつけていなかったのですが
昨日の授業で友人と悩んでいて結局つけていました。
範囲については上でご指導いただいたとおりにしたいと思っています。
本当にありがとうございました。

お礼日時：2007/12/01 01:29