シンプソン則について

大学の課題で　∫0～5　e^(-3t) * cos^2 *(2t) をシンプソン則を用いて求めるプログラムを作成せよ。積分の刻みは0.01毎　　という課題がでました。
私は授業で作ったプログラムを元に、以下のプログラムを作ったのですが、答えが大きすぎて（5より大きいと間違いと言っていました）受け取ってもらえません。どこがいけないのでしょうか？
あと、ｎには　どのような数字をいれればよいのでしょうか？ｎは使わないのでしょうか。
お願いします。

#include "stdio.h"
#include "math.h"

void main()
{
int n;
double i, a, b, h, S, S1, S2, S3;
double f (double);
a = 0.;
b = 5.;

printf ("N=");
scanf ("%d", &n);
h = (b - a) / (double)n;
S1 = 0.;
S2 = 0.;
S3 = f (a) + f (b); /* 両端の値*/
/* 奇数番目の値*/
for (i = 0.01; i < 5; i = i + 0.02)
{
S1 = S1 + 4. * (f (a + (double)i * h));
}
/* 偶数番目の値*/
for (i = 0.0; i < 5; i = i + 0.02)
{
S2 = S2 + 2. * (f (a + (double)i * h));
}
S = (h / 3.) * (S1 + S2 + S3);
printf ("%d %lf\n", n, S);
}
double f (double t)
{
double y;

y = (exp(-3. * t) / 2.) + ((exp(-3.* t) * cos(4. * t)) / 2.);
return y;
}

No.10 回答者： aquarius_hiro 回答日時： 2007/07/12 22:29

A No.8 です。

アルゴリズムはそれで合っていると思いますよ。
（細かいところは見ていないですが…。）

ただ i はもともと double で宣言しているので、(double)i としなくて良いと思います。そのまま a + i で良いです。

No.9 回答者： asuncion 回答日時： 2007/07/06 22:51

こうだと思います。

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

double f(double x)
{
　　return exp(-3 * x) * cos(2 * x) * cos(2 * x);
}

int main(void)
{
　　double a, b, h, s;
　　int n, i;
　　
　　a = 0, b = 5;
　　h = 0.01;
　　n = (b - a) / h + 0.5;　　　　// +0.5 は丸め誤差への対応
　　for (s = 0, i = 1; i < n; i++)
　　　　s += ((i % 2 == 0) ? 2 : 4) * f(a + i * h);
　　s += f(a) + f(b);
　　s *= h / 3;
　　printf("面積=%f\n", s);
　　return 0;
}

（注）インデントのため、全角空白を使っています。

No.8 回答者： aquarius_hiro 回答日時： 2007/07/06 22:08

だんだんまどろっこしくなってきたので…

もう言ってしまいますね。。(^^;)

i は 0.02 ずつ増えています。

a + (double)i*h は、例えば偶数のほうでは、

a, a + 0.02*h, a + 0.04*h, .....

と増えていきますね。

しかし、ほんとうは、a, a + 0.02, a + 0.04, ... でよいのですから、*h が余分ですよね。

a + i だけで良いです。

奇数のほうでも同様です。


もし a + (double)i * h の形にしたければ、

h = 0.01 刻みで増えていくようにするには、i は 2 刻みで増やさないといけませんね。（i=1,3,5,... と i=0,2,4,... ）

この回答へのお礼

なかなか理解出来ずにすみません。

このプログラムでよいのでしょうか？
本当にありがとうございました！
#include "stdio.h"
#include "math.h"

void main()
{
int n;
double i, a, b, h, S, S1, S2, S3;
double f (double);
a = 0.;
b = 5.;

printf ("N=");
scanf ("%d", &n);
h = (b - a) / (double)n;
S1 = 0.;
S2 = 0.;
S3 = f (a) + f (b); /* 両端の値*/
/* 奇数番目の値*/
for (i = 0.01; i < 5; i = i + 0.02)
{
S1 = S1 + 4. * (f (a + (double)i));
}
/* 偶数番目の値*/
for (i = 0.0; i < 5; i = i + 0.02)
{
S2 = S2 + 2. * (f (a + (double)i));
}
S = (h / 3.) * (S1 + S2 + S3);
printf ("%d %lf\n", n, S);
}
double f (double t)
{
double y;

y = (exp(-3. * t) / 2.) + ((exp(-3.* t) * cos(4. * t)) / 2.);
return y;
}

お礼日時：2007/07/06 22:18

No.7 回答者： aquarius_hiro 回答日時： 2007/07/06 21:51

そうですね、i は 0.02 ずつ増えるわけですね。

そのとき、

f() の引数

a + (double)i * h

は、a = 0 （の近く）から、どんなふうに増えていくと思いますか？

具体的に数字を順に書いてみたらどうでしょうね？

b = 5 （の近く）まで変化していかないといけないわけですが…。

この回答へのお礼

a + (double)i * h の aも0.02ずつ値が大きくなるということでしょうか？つまり、　i + (double)i * h でしょうか？

お礼日時：2007/07/06 21:57

No.6 回答者： asuncion 回答日時： 2007/07/06 21:03

> 半角の公式や余弦定理などを使い　分解してみました！

シンプソンの公式にはそういう変形が必要だとは書いていないです。

No.5 回答者： asuncion 回答日時： 2007/07/06 20:29

> ∫0～5　e^(-3t) * cos^2 *(2t)

積分したい関数と

> y = (exp(-3. * t) / 2.) + ((exp(-3.* t) * cos(4. * t)) / 2.);

実際のコードに書かれたこの式との関係がわかりません。

この回答へのお礼

半角の公式や余弦定理などを使い　分解してみました！

お礼日時：2007/07/06 20:42

No.4 回答者： aquarius_hiro 回答日時： 2007/07/06 20:22

できるだけ自分で考えてほしいので、まどろっこしい書き方になってしまいますが…。

> double fの中での　t にうまく代入されていないということですか？

そういうことではなく…

for (i = 0.01; i < 5; i = i + 0.02)

では、i はどのように増えていっていますか？

そして、

f(a + (double)i * h) と書かれていますが、

この変数は意図したように値が増えていっていますか？

よ～く、考えてみてください。

h = ?

もしこの式なら、i はどう増えていくべきですか？

修正の仕方は二通りあります。

この回答へのお礼

for (i = 0.01; i < 5; i = i + 0.02)は奇数だけを考えるようにしたので0.02ずつ増えていくようにしました。その代わり偶数を考えるfor文も作り、こぼれがないようにしたつもりです。

f(a + (double)i * h) の式の意味が分からなくなってしまいました。今回の問題には使えないような気がしてきました・・・。

お礼日時：2007/07/06 20:38

No.3 回答者： Tacosan 回答日時： 2007/07/06 19:46

あ, h はその意味なんですね. それなら 3 ですね. 了解です.

あと, ちょっとそのプログラムを動かしてみたんだけど危険なところがいくつかありますね. 簡単に列挙すると
・stdio.h や math.h は #include <～> でインクルードすべし
・f の宣言は main の中でしない方が安全
・やたらと無駄なキャストをし過ぎ
・for で double の変数を動かすのは危険
くらい?
ついでにいうと, #2 で挙げられたコメントとは別件で, 多分もう 1つ間違ってますね. 「偶数番目の値」のコメントのあとの for 文があやしい感じ.

この回答へのお礼

y = (exp(-3. * t) / 2.) + ((exp(-3.* t) * cos(4. * t)) / 2.)　のtに値を入れていくのはdouble fのなかで　forで回せばいいのですか？

for (i = 0.02; i < 5; i = i + 0.02)にすればよいでしょうか？
質問ばかりですみません！

お礼日時：2007/07/06 20:05

No.2 回答者： aquarius_hiro 回答日時： 2007/07/06 19:13

どこが間違っているのか指摘したら、getherさんの勉強にならないと思うので、ヒントをお教えします。

きざみ幅が0.01で、区間の幅が5なら、きざみの数 n はどうなるか考えてみてください。また、h = (b-a)/n と書いた h は何を表しているのかも考えてみてください。0.01とhの両方がプログラムの中に書かれていますが…。

そして、f(x) の x は loop の中で思ったように増えていっているか、よく考えてみてください。

この回答へのお礼

ｎ＝500ということでしょうか？
ｈは0.01刻みごとの幅だと考えています。

double fの中での　t にうまく代入されていないということですか？
回答ありがとうございました！！

お礼日時：2007/07/06 19:23

No.1 回答者： Tacosan 回答日時： 2007/07/06 18:55

あれ? シンプソン則って, 3 で割るんだっけ? なんか, 6 だったような気がするんだけど....

この回答へのお礼

ネットで調べたところ　3で正しいようです。
迅速な回答ありがとうございます。

お礼日時：2007/07/06 19:00