fortranで図形の面積を求めるには？

今fortranのプログラムをやっているのですが、
図形の面積を求めたいのですがどうやったらいいのでしょうか？
とりあえず小さな区間に区切って足し算を行うようにしているのですが、うまく動きません・・・。
ＤＯループをまわし過ぎなのでしょうか？？？
一応下にプログラムの一部を書いておきます。

DO 444 J= 1,5000
T=FLOAT(J)*0.1

IF (T.GE.281)THEN
GL3 = -0.01105*T+4.104972
ELSE
　　 GL3 = -0.01105*T-2.1050
ENDIF
　　　　・
　　　　・
　　　　・
444 CONTINUE

誰か助けてください！
宜しくお願いします。

No.6 ベストアンサー 回答者： a-kuma 回答日時： 2001/07/04 10:02

質問者本人抜きの議論は、ここでは御法度ですが、質問者の理解の助けに

なると信じ、実数の比較の話に決着をつけます。

ubonoti01 さんの理解は、理由を抜きに「実数の比較が危ないらしい」と
いう風になっているのがまずい。

今回の質問のプログラム

DO 444 J= 1,5000
　T=FLOAT(J)*0.1
　IF (T.GE.281)THEN

であれば、EQUAL を含めて、問題なく判定されます。問題が出るのは、

T = 0.0
DO 444 J = 1, 5000
　T = T + 0.1
　IF (T .GE. 281.0) THEN

となっている場合。でも、同様の書き方でも

T = 0.0
DO 444 J = 1, 5000
　T = T + 0.5
　IF (T .GE. 281.0) THEN

であれば、問題なく動く計算機が *多い* です。

どうしてこんなことになるかというと、（Fortran で扱う）実数が二進数で
表現されている（正確にいえば、二進数の指数表現）からです。

0.1 は、十進数で書くと切りが良いですが二進数で表現すると無限小数に
なります。だから二進数の世界では

　　0.1 × 10 ≠ 0.1 ＋ … ＋ 0.1 (10回足す)

となってしまいます。

先ほど例に出した 0.5 は、二進数の世界でも無限小数ではありません。
だから、

　　0.5 × 10 ＝ 0.5 ＋ … ＋ 0.5 (10回足す)

が成立します。

で、このへんの話は「実数を二進数で表現すること」の問題なので、
Fortran 固有の問題ではなく、C でも同様の問題を抱えます。

bob> 実数の大小比較に信頼性がない

のではなく、演算結果が等しいかどうかを判断するときに実数演算の
誤差を考慮しなくてはいけない、と言うことです。

ubonoti01> IF (J.GE.2810)THEN の方が好ましい

とあるのは、その「考慮」のやりかたのひとつです。



(補足 その１)
先程、0.5 を足して行くケースで「問題なく動く計算機が *多い* です」と
いう表現を使いました。それは、実数の表現が二進数表現ではない計算機が
あるからです。いわゆる IBM Format と呼ばれる実数表現では 16進数の
指数形式で表現されます。

なので、多分 0.5 では 0.1 と同様の問題が出るはずで、0.0625 (1÷16)なら
問題は発生しません。

(補足 その２）
COBOL では、ふつう実数は BCD コードで桁数を指定して扱います。なので、
オーバフローやアンダーフローが無いように桁を決めてあげれば、こういった
問題は出ません。

No.5 回答者： bob 回答日時： 2001/07/04 00:27

便乗質問になってしまいますが、私は基本的にCプログラマでFortranでの計算はあまりやったことがないのでFortran独特のくせは確かによく知らないので、一応質問させていただきます。

実数の大小比較に信頼性がないということは、そのままとるとFortranは実数が使えないというのと同義です。主に科学技術計算に使われてきたFortranの歴史からいくとちょっと信じがたいのですが、これは本当なのでしょうか？

No.4 回答者： ubonoti01 回答日時： 2001/07/03 23:30

再度ubonoti01です。

わたしがIF (T.GE.281)THEN は、IF (J.GE.2810)THEN の方が好ましいと申したのは、
・IF文の中でREAL同士を比較するのは危険と言う意味です。IF文の中での値の比較はINTEGER同士の比較であるべきです。出所が別のREAL同士の比較では、本来一致する筈が真でなく偽と判定されるケースが多いからです。これはコンピュータ内部での数値の扱いに起因します。

No.3 回答者： bob 回答日時： 2001/07/03 22:54

やはり、最初に問題の説明が不足していることを指摘しておきます。

この質問内容では推測による回答しかできません。

で、推測に基づく回答ですが、やっていることは数値積分ですね。おそらく１次式で表わされる２直線とx軸（t軸）とy軸（gl3軸）で囲まれる部分の面積を求めようとしているのだと思われます。

で、予想される問題点ですが、多分変数Tに掛けている定数-0.1105のうちどちらか（多分後者）は符号が逆でしょう。切片も-2.1050ではなくて2.1050では？

if分岐の条件についてですが、ubonoti01さんと意見が割れますが、この場合は積分変数による分岐なのでTで分岐する方で良いと思われます。
おそらく意味的にもっといいのは各式を関数化しておいて、IF([２式].GE.[１式])THENと書いてしまうことです。計算量は若干多くなりますが、281という一見意味不明なパラメータを直接書く必要はなくなります。

そして、おそらく一番決定的な問題は積分の和をとってないことです。先頭にS=0をいれて、CONTINUEの前にS=S+GL3を入れるか、GL3を初期化しておいて２つの式をGL3=GL3+...にする必要があると思われます。

最後に、ループの終了条件がJ=5000というのはちょっと怪しいと思われます。この時GL3はゼロを割り込んでいます。おそらく3715以下の数値になると思われます。

No.2 回答者： ubonoti01 回答日時： 2001/07/03 14:55

もう少し背景なり状況なりを示した質問をすべきと思いますが・・・・・・。

プログラムだけを見た範囲では、a-kumaさんのご指摘のほかに、
IF (T.GE.281)THEN
は、
IF (J.GE.2810)THEN
の方が好ましいと思います。

No.1 回答者： a-kuma 回答日時： 2001/07/03 14:21

どう上手くいかないんですかね？

コンパイルができない、期待している数字と結果が違う、計算に時間がかかり
過ぎる、etc ．．．

このコードの断片だけを見ると、コンパイルはできそうだし、ぎざぎざの
のこぎりのような形の面積を求めたいのかな、位しか想像つきません。

もしそうだとしたら、GL3 の初期化がされていないところと、Ｔの判断条件
のところが怪しいかな、くらいを思い付きます。

　　　GL3 = 0.0 ← ★ 必要なし？
　　　DO 444 J= 1,5000
　　　　T=FLOAT(J)*0.1
　　　　
　　　　IF (T.GE.281)THEN ← ★ 281？ 28.1？ 不等号の向きは合ってる？
　　　　　GL3 = -0.01105*T+4.104972 ← ★ 右下がりで
　　　　ELSE
　　　　　GL3 = -0.01105*T-2.1050 ← ★ こっちも右下がり？
　　　　ENDIF ← ★ 終わっちゃっていいの？
　　　　　　　　・ ← ★ 続いているみたいだけど
　　　　　　　　・
　　　　　　　　・
　444 CONTINUE

ぎざぎざな形の面積だったら、こんな感じになるんじゃない？

IF (T .LT. 100.0) THEN
　GL3 = ...
ELSE IF (T .LT. 200.0) THEN
　GL3 = ...
ELSE IF (T .LT. 300.0) THEN
　・
　・
ENDIF