大学で現在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 ***
＊＊＊ここまで＊＊＊
もしお分かりになる方がいらっしゃいましたら、
お手数ですがご指導いただけないでしょうか？
お願いいたします。
大学で現在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 ***
＊＊＊ここまで＊＊＊
もしお分かりになる方がいらっしゃいましたら、
お手数ですがご指導いただけないでしょうか？
お願いいたします。
モンテカルロ法で円周率の推定値を計算することを最近習ったのですが、定積分でもそれが可能なのを知り、どうやってプログラムを組めばいいのか分からず、困っています。
例えば、定積分∫[0→1]x^2dx=1/3～0.333（[0→1]というのは、積分範囲です。）をモンテカルロ法で計算すると、どういうプログラムを組めばいいのでしょうか？
わかる範囲で書いてみたのですが…積分の範囲をどうやってプログラミングすればいいのか、いまいち分かりませんでした。
教えていただけると、助かります。よろしくお願いします。
RANDOMIZE
INPUT n
SET WINDOW -0.1,1.1, -0.1,1.1
DRAW GRID
SET POINT STYLE 1
LET sumin=0
FOR i=0 TO n
LET x=RND
LET y=RND
SET POINT STYLE 2
IF y<x*x THEN
SET POINT COLOR 4
LET sumin=sumin+1
END IF
! PRINT USING "(%.####, %.####)": x,y
PLOT POINTS: x,y
NEXT i
PRINT 1*1*sumin/n
END
モンテカルロ法で円周率の推定値を計算することを最近習ったのですが、定積分でもそれが可能なのを知り、どうやってプログラムを組めばいいのか分からず、困っています。
例えば、定積分∫[0→1]x^2dx=1/3～0.333（[0→1]というのは、積分範囲です。）をモンテカルロ法で計算すると、どういうプログラムを組めばいいのでしょうか？
わかる範囲で書いてみたのですが…積分の範囲をどうやってプログラミングすればいいのか、いまいち分かりませんでした。
教えていただけると、助かります。よろしくお願いします。
RANDOMIZE
INPUT n
SET WINDOW -0.1,1.1, -0.1,1.1
DRAW GRID
SET POINT STYLE 1
LET sumin=0
FOR i=0 TO n
LET x=RND
LET y=RND
SET POINT STYLE 2
IF y<x*x THEN
SET POINT COLOR 4
LET sumin=sumin+1
END IF
! PRINT USING "(%.####, %.####)": x,y
PLOT POINTS: x,y
NEXT i
PRINT 1*1*sumin/n
END
