C言語 無限区間積分の数値計算のプログラムのエラー

Question

###前提・実現したいこと
次のような課題に取り組んでいます。

無限区間積分 ∫_0^∞ (exp(-x)) dx をスケール変換 x=(1+u)/(1-u) を用いて,有限区間の定積分に変換し、積分値を求めよ。

スケール変換を行った後の積分の式は
∫_(-1)^(1) (exp{-(1+u)/(1-u)}*{2/(1-u)^2}) du
になると考えています。

シンプソンの公式を用いたプログラムを組んでみて結果を出力しようとすると、-1.#IND00 と出てきてしまい、ちゃんと値が計算されていません。おそらく原因は、どこかで0で割る計算が行われてしまっていることだと推測しているのですが、どこが原因なのかがわかりません。

ソースコードの間違いを指摘していただけると大変助かります。
よろしくお願いいたします。

###発生している問題・エラーメッセージ

コンパイルエラーではないのですが、
出てくる数値が以下のようになってしまいます
-1.#IND00

###該当のソースコード

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double f(double x){
	double y;
	y = (1/exp((1+x)/(1-x)))*(2/((1-x)*(1-x)));
	return y;
}

int main(void){
	int n;
	double h, I, sum1 = 0, sum2 = 0;

	h = (1 - (-1)) / 1000;

	for (n = 1; n < 1000 / 2 - 1; n++){
		sum1 += f(-1 + 2 * n*h);
	}

	for (n = 1; n < 1000 / 2; n++){
		sum2 += f(-1 + (2 * n - 1)*h);
	}

	I = (h / 3)*(f(-1) + f(1-(2/1000)) + 2 * sum1 + 4 * sum2);

	printf("%f\n",I);
	
	return 0;
}

Answer 1

どこに問題があるかを突き止めるには、まずはデバッグ文です。
以下のようにしデバッグ文を入れてみました。

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double f(double x){
    double y;
    y = (1/exp((1+x)/(1-x)))*(2/((1-x)*(1-x)));
    printf("x=%lf,y=%lf\n",x,y);
    return y;
}

int main(void){
    int n;
    double h, I, sum1 = 0, sum2 = 0;

    h = (1 - (-1)) / 1000;
    printf("h=%f\n",h);

    for (n = 1; n < 1000 / 2 - 1; n++){
        sum1 += f(-1 + 2 * n*h);
    }
    printf("sum1=%f\n",sum1);

    for (n = 1; n < 1000 / 2; n++){
        sum2 += f(-1 + (2 * n - 1)*h);
    }
    printf("sum2=%f\n",sum2);

    I = (h / 3)*(f(-1) + f(1-(2/1000)) + 2 * sum1 + 4 * sum2);

    printf("I=%f\n",I);

    return 0;
}

これでそれぞれの計算結果、関数の引数などがどのような値になっているかわかります。
結果がこれです。

h=0.000000
x=-1.000000,y=0.500000
x=-1.000000,y=0.500000
～省略～
x=-1.000000,y=0.500000
x=-1.000000,y=0.500000
sum2=249.500000
x=-1.000000,y=0.500000
x=1.000000,y=-nan
I=-nan

問題は3箇所。

• hが0になっている
• 関数fに同じ値、戻り値が同じ値
• 一番最後の関数fの戻り値がnan

いずれも理由は同じです。
以下の計算式が整数同士の計算なので、小数点以下がまるめられてしまっているからです。

h = (1 - (-1)) / 1000;
I = (h / 3)*(f(-1) + f(1-(2/1000)) + 2 * sum1 + 4 * sum2);

これをキャストなり小数値を明示的に使うことで、まるめないようにすることで対処できます。

h = (1 - (-1)) / 1000.0;
I = (h / 3)*(f(-1) + f(1-(2/1000.0)) + 2 * sum1 + 4 * sum2);

で、結果がこれ。

h=0.002000
x=-0.996000,y=0.501001
x=-0.992000,y=0.502004
～省略～
x=0.990000,y=0.000000
x=0.994000,y=0.000000
sum2=250.000042
x=-1.000000,y=0.500000
x=0.998000,y=0.000000
I=1.000000

とりあえずnanは免れましたが、1というのは想定通りでしょうか。
計算式についての証明は行っていませんので、あとは自力解決をお願いします。

Comments

[jackie687456]

まず、デバッグというもののやり方を学ばせていただいてありがとうございます。できるだけ自分の力で見つけられるように努力していきます。

質問に対する回答についてですが、非常にわかりやすい説明で大変助かりました。整数同士の計算式だと、小数点以下が無視されてしまうんですね…。そんなところにミスがあったとは。質問欄に誤った推測を載せてしまい申し訳ありませんでした。

無事、こちらでも数値を出すことができました。１で想定通りです。
本当にありがとうございました。