Doループとファンクション（fortran)

質問

課題で初めてプログラムに取り組んでいるものです。よろしくお願いします。

台形則を使って楕円の面積（の４分の１）を求めるプログラムを作ったのですが、ファンクションで関数を指定するとうまくいきません。

具体的には
func(x)=4.0*sqrt(1.0-x**2.0)
という関数で、0<=x<=1を１０００万分割して面積を近似するという課題です。
本来ならfunc(x)は徐々に減少し、x=1で０になるはずが、実際に走らせてみると、それより手前で０になってあとはNaNがでます。

func(x)以外（x座標の分割の仕方等）はうまく行っており（確認済み）、function形式をやめてメインのプログラムに直接式を書き込んだらうまくいきました。

自分ではなんでこんなことが起こるのか皆目見当が付きません。ご教授お願いしますm(_ _)m

回答No.2

koma1000nin · ベストアンサー

a，bそれぞれを不動小数点数とすると、

a^b = e^(b*log(a))

として計算されます。底はeです。
この計算だけでも、EXP(X)の計算，ALOG(X)の計算，乗算を経由するたびに誤差が積もリます。そのうえ、単精度浮動小数点数と倍精度浮動小数点数の精度の違いで、誤差がいっそう積もりますので、べき計算は"X**2"とすべきでしょう。
関数で計算する場合は、ご自身で倍精度計算に変更してみるのはいかがでしょう。

4.0*sqrt(1.0-x**2.0)
⇒DBLE(4.0)*DSQRT(DBLE(1.0)-X**2)

引数Xは倍精度とします。関数FUNNCも倍精度宣言すれば、メイン側で計算するのと同じ結果が得られるでしょう。

お礼

整数以外の結果が出る計算をするときは全ての数を浮動小数点にしなければならないと思っていました(^^;;
誤差というものをもっと学ぶ必要がありますね・・・。

数値の倍精度指定の仕方まで教えてくださって有難うございました!!

回答No.1

ultraCS

課題だと答えを出すのは為にならないので

たぶん、誤差の問題ですね。
何も宣言しないと、XやFUNC、1.0、2.0は単精度になっているはずです。単精度だと、有効数字の範囲が、1000万をフルに表現できるだけありません。
ではどうするかは、ちょっと考えてみてください。

数値計算に取り組むなら、まず、誤差について勉強してみましょう。

なお、余談ですが、べき乗は処理系にもよりますが、X**2.0とX**2、X*Xでは内部の処理が違います、べき数を浮動小数点にすると、一度対数を呼ぶので計算課程が複雑になり、精度、速度が下がることがあります。2乗であれば、x**2かX*Xと書きましょう。

お礼

有難うございました!!
数値を倍精度にする方法はまだ知らないので、自分で探してみます(^^)