マクローリン展開による三角関数の近似値計算プログラム

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
アバター
kerotan0820
記事: 91
登録日時: 13年前
住所: 東京都
連絡を取る:

マクローリン展開による三角関数の近似値計算プログラム

#1

投稿記事 by kerotan0820 » 10年前

今現在、課題でマクローリン展開による三角関数の近似値計算をC言語で書いております。

仕組みとプログラムの構成自体も概ね理解し組むことが出来たのですが答えが微妙に合いません。

コード:

/************************************/
/*  プログラミング演習レポート#00  */
/* 【 ここにレポートの題目を書く 】 */
/*                  学籍番号: */
/*                     作成者:*/
/*               作成日:2014/05/20 */
/************************************/

/************************************************
		ヘッダー、定義・マクロ
*************************************************/
#include <stdio.h>
#include <math.h>

#define PI 3.141592653589793


/************************************************
			プロトタイプ宣言
*************************************************/
int menu(void);
void msin(void);
void mcos(void);
void mexp(void);
double macsin(double rad, int n);
double maccos(double rad, int n);
double macexp(double x, int n);

/************************************************
				メイン
*************************************************/
int main( ) {
	int num;
//	setbuf(stdout,0);//msPEwb2Kの制限:main関数の最初の実行文の前に置く
	do{
		switch( num = menu() ){
			case 1:
				msin();
				break;
			case 2:
				mcos();
				break;
			case 3:
				mexp();
				break;
			case 4:
				printf("終了します\n");
				break;
			default:
				printf("1~4の半角数字を入力してください\n");
		}
	}while(num != 4);
}

/************************************************
				メニュー
*************************************************/
int menu(void){
	int n;
	printf("-----------------------\n"
		   " 1:sine 関数        \n"
		   " 2:cosine 関数      \n"
		   " 3:exponetial 関数  \n"
		   " 4:終了             \n"
		   "-----------------------\n");
	scanf("%d", &n);
	return n;
}

/************************************************
				sin関数 出力
*************************************************/
void msin(void){
	int x;
	double M, m5, m6, m7, e5, e6, e7, rad;

	puts("-------------------------------------------------------------------------------------------");
	puts("角度 | 標準関数の値| 5項まで   相対誤差|  6項まで   相対誤差| 7項まで    相対誤差");
	puts("-------------------------------------------------------------------------------------------");
	for(x = 0; x <= 180; x += 5){
		rad = x * PI / 180.0;
		M = sin(rad);
		m5 = macsin(rad, 5);
		m6 = macsin(rad, 6);
		m7 = macsin(rad, 7);
		if(fabs(M) > 0.01){
			e5 = fabs( (M - m5) / M );
			e6 = fabs( (M - m6) / M );
			e7 = fabs( (M - m7) / M );
		}else{
			e5 = e6 = e7 = 0.0;
		}
		printf("%5d|%13.10f|%13.10f %9.1e|%13.10f %9.1e|%13.10f %9.1e\n", x, M, m5, e5, m6, e6, m7, e7);
	}
}

/************************************************
				sin関数 計算
*************************************************/
double macsin(double rad, int n){
	int i, sign, bunbo;
	double ans, tmp;

	tmp = rad;
	ans = rad;
	sign = 1;
	bunbo = 1;

	for(i = 1; i <= n; i++){
		sign *= -1;						//符号反転
		tmp *= rad * rad;				//奇数乗(3 -> 5 -> 7..)
		bunbo *= ( 2 * i + 1 ) * ( 2 * i );
		ans = ans + sign*(tmp / (double)bunbo);
	}
	return ans; 
}
(少しプログラムが長かったため不必要な部分は消しました)
次の画像が出力例として提示されたものと、私の実行結果です。

非常に近い数字が出ているのでいったい何が原因となっているのか分からず困っております。
repo.png
CCF20140520.jpg
ご教示いただけると幸いです。
けろけろにゃー (」・ω・)」うー!

かずま

Re: マクローリン展開による三角関数の近似値計算プログラム

#2

投稿記事 by かずま » 10年前

kerotan0820 さんが書きました:非常に近い数字が出ているのでいったい何が原因となっているのか分からず困っております。
マクローリン展開は、x = 0 の近傍での近似値なので
x が 0 から遠ければ遠いほど誤差が大きくなります。

for(x = 0; x <= 180; x += 5){ を
for(x = 360; x <= 360+180; x += 5){ に変更してみてください。
もっと誤差が大きくなることがわかるでしょう。

旅路のきのこ

Re: マクローリン展開による三角関数の近似値計算プログラム

#3

投稿記事 by 旅路のきのこ » 10年前

計算している項の数が一つ多いかもしれません。
出力例の6項までの値と、kerotan0820さんの5項までの値がおんなじです。

あと、きっちり確かめてはいませんが、ループ内のiの値が大きくなると、bunboの値がおかしくなってしまっています。
32ビットに入りきらないんでしょうか?
macsin(rad, 10);とかやってみましたが、凄い値が出てきます。
誤差がどんどん大きくなっていってるんですね、多分。
参考になれば幸いです。

アバター
みけCAT
記事: 6734
登録日時: 13年前
住所: 千葉県
連絡を取る:

Re: マクローリン展開による三角関数の近似値計算プログラム

#4

投稿記事 by みけCAT » 10年前

細かい実数の計算は、使用するCPUや計算順序、最適化などの条件により、誤差が出る(違った値になる)可能性があります。
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

watcher

Re: マクローリン展開による三角関数の近似値計算プログラム

#5

投稿記事 by watcher » 10年前

 最初は計算機誤差なのかなと思いましたが,角度が小さいと高次の項が十分小さいので見えてこないだけで,計算機誤差では説明できない明らかな違いがありますね.
旅路のきのこさんの指摘されたように,「5項まで」 と書いてある列は第 6 項 (x の 11 乗の項) まで計算した結果ですね.
そうなると,出力例と kerotan0820 さんの結果は列をずらせば同じになるはず,と思いましたが,
第 7 項 (x の 13 乗) 以上の計算では別の問題が起きているようですね.
 関連するこんな質問がありました.
http://detail.chiebukuro.yahoo.co.jp/qa ... 1320627493

 私も勉強になりました.

アバター
GRAM
記事: 164
登録日時: 13年前
住所: 大阪

Re: マクローリン展開による三角関数の近似値計算プログラム

#6

投稿記事 by GRAM » 9年前

オーバーフローに関してはほかの人が指摘しているとおりでしょう。
というかこの手の計算でint型を使う必要は皆無だと思います。
大体多くのdouble型は仮数部が50bit以上あります。すなわちdouble型はint型の扱う範囲であれば誤差なく表現できます。
多倍長整数が必要になるような大きな整数を扱うのでない限り、doubleを使えばいいと思います。
このことにつけ加えるなら誤差の原因は

①おそらく標準の実装はforループを使わないで、誤差が最小になるようにあらかじめ計算順序を最適化しているから
この辺は計算機数学の話で、どういう計算をすれば誤差が小さくなるのかはsinを求めるアルゴリズムとは本質的には関係ない

②sinもcosも45°の値までを用いればよい。というか基本的に0近傍以外の値はあまり信用できない。
sin50°= cos40° 表からも明らかなように45°までを用いれば誤差は10^-10 以下のオーダーになる。
そしてほとんどの場合において、実用上10^-10オーダーというのは話にならないくらい小さい誤差といえる。

アバター
kerotan0820
記事: 91
登録日時: 13年前
住所: 東京都
連絡を取る:

Re: マクローリン展開による三角関数の近似値計算プログラム

#7

投稿記事 by kerotan0820 » 9年前

たくさんの解答ありがとうございます。

かずまさま>>
確かに誤差が大きくなりました。
提出時にはそのようなコメントを添えておこうと思います。
ありがとうございました。


旅路のきのこさま>>
まったく気がつきませんでした
はじめから ans に rad を入れてスタートしていたので確かに1つ項を多く計算しておりました。

あと、オーバーフローはその通りのようでした。
一応64bitPCではありますが・・・・。
分母での15の階乗が既に1兆?を上回る数で、それが原因のように感じましたがマクローリン展開をプログラムにするという課題はクリアできましたのでこの状態で一度提出したいと思います。
非常に勉強になりました。
ありがとうございました。


みけCATさま>>
非常に興味深いです。
今まで深く変数のサイズについて気にすることが無かったので非常に良い勉強になりました。
ありがとうございました。

watcherさま>>
リンクまで添えていただきありがとうございます。
確かに結果がずれていました。
修正したところ、7項の問題以外はクリアできました。
興味深い話だったので、課題提出後調べてみようと思います。
ありがとうございました。


GRAMさま>>
詳しくありがとうございます。
intからdoubleに分母の階乗部分を計算したら誤差が解答例と同じ範囲に収まりました!

確かに 45°までの値を用いれば良いですね。どの手を打つべきか先生に尋ねてみたいと思います。
ありがとうございました。


皆様ありがとうございました。
けろけろにゃー (」・ω・)」うー!

かずま

Re: マクローリン展開による三角関数の近似値計算プログラム

#8

投稿記事 by かずま » 9年前

kerotan0820 さんが書きました: 分母での15の階乗が既に1兆?を上回る数で、それが原因のように感じましたがマクローリン展開をプログラムにするという課題はクリアできましたのでこの状態で一度提出したいと思います。
オーバフローして間違っていることが分かっているのに、それで提出ですか?
double bunbo; にするだけじゃないですか。
最初の私の回答が的外れだったことをお詫びします。

アバター
kerotan0820
記事: 91
登録日時: 13年前
住所: 東京都
連絡を取る:

Re: マクローリン展開による三角関数の近似値計算プログラム

#9

投稿記事 by kerotan0820 » 9年前

かずま さんが書きました:
kerotan0820 さんが書きました: 分母での15の階乗が既に1兆?を上回る数で、それが原因のように感じましたがマクローリン展開をプログラムにするという課題はクリアできましたのでこの状態で一度提出したいと思います。
オーバフローして間違っていることが分かっているのに、それで提出ですか?
double bunbo; にするだけじゃないですか。
最初の私の回答が的外れだったことをお詫びします。

すみません、かずまさんへの返信にdouble型にしたことで解決したことを記述できておりませんでした。
これからこのような大きな値を用いるプログラムが増えそうなので気をつけます。
ありがとうございました。
けろけろにゃー (」・ω・)」うー!

閉鎖

“C言語何でも質問掲示板” へ戻る