sin(x)の近似について

sinｘの数値計算
任意のｘに対するsinxの値をマクローリン展開を利用して近似し、誤差の限界（ｎ番目の値が1*10^-8）になるまでもとめよ。


という問題なんですが、for文でいろいろやってみたのですが、ｎ番目の値が1*10^-8までというのがどうしてもできません。
C言語です。
たぶん、階乗の部分が間違ってると思われます。
下のが自分で考えたものです。アドバイスいただけたら嬉しいです。



#include <stdio.h>

int main(void)
{
signed int i,v;
double x=-1;
double a,m,n=1.0;
double sinx = 0.0;


a=x;
v=1;
printf("a=%f\n",a);
printf("sinx = %f\n",sinx);

while(1)
{
sinx = sinx + a;
m=x;
m=m*x;
a = (-1)*v*m/((n+1)*(n+2));
n=n+2;
v=v*(-1);
if ((a < 1e-8) && (a > -1e-8)) break;



}
printf("sinx = %f",sinx);


return (0);

}

No.2 ベストアンサー 回答者： chie65535 回答日時： 2009/07/16 11:30

#include <stdio.h>

/*#include <math.h>*/ /*検算用に使うsin(x)用。検算が必要ならコメントを外す*/

int main(void)
{
int n,v;
double a,k,l,xx;
double x;
double sinx;

printf("sin(x)のxを入力:");
scanf("%lf",&x);

/*誤差の限界は「10の-8乗」なので精度を「8桁」に指定*/
/*printf("sin = %.8f\n",sin(x));*/ /*検算用にライブラリ関数でのsin(x)を表示。検算が必要ならコメントを外す*/

/*変数の初期化*/
sinx = 0.0;
n = 1; /*２ｎ＋１。1,3,5,7…と増加*/
v = 1; /*加算か減算かのフラグ。1,0,1,0,1,0…と変化*/
k = 1.0; /*（２ｎ＋１）！。最初は１！なので１*/
l = x; /*ｘの（２ｎ＋１）乗。最初はｘの１乗なのでｘ*/
xx = x * x; /*ｘの１乗にｘの２乗を掛け、ｘの３乗、５乗、７乗…を作る為に用意した「ｘの２乗」*/
while(1)
{
a = (1.0 / k) * l; /*赤枠の部分式*/
if (v) sinx += a; /*ｎが1,5,9…なら加算*/
else sinx -= a; /*ｎが3,7,11…なら減算*/
n += 2; /*nを2づつ加算*/
k *= (n - 1) * n; /*（２ｎ＋１）！を計算*/
l *= xx; /*ｘの（２ｎ＋１）乗を計算*/
v = !v; /*加算と減算の入れ替え*/
/*誤差の限界判定はsinxにaを加減算した後で行う*/
/*そうしないとxが3.1415926535の時に結果が*/
/*ライブラリのsin関数と食い違ってしまう*/
if ((a > -1e-8) && (a < 1e-8)) break; /*誤差の限界*/
}
/*誤差の限界は「10の-8乗」なので精度を「8桁」に指定*/
printf("sinx = %.8f",sinx); /*結果を表示*/
return (0);
}

この回答へのお礼

たびたび、質問に答えていただきありがとうございました。感謝しています。

お礼日時：2009/07/16 14:41

No.4 回答者： hitokotonusi 回答日時： 2009/07/16 12:34

もう一度失礼します。

少しの間違いを除けば基本的に質問者さんのコードであっているし、
そちらのほうがセンスがいいので、それをもとにしたほうがベターです。

質問者さんはわかっているはずなので念のためですが、サインの展開式

sin x = Σ[m=0,∞] (-1)^m x^{2m+1}/(2m+1) ! = Σ[m=0,∞] a(m)

から

a(m) = (-1)^m x^{2m+1}/(2m+1) !
　　= (-1)^(m-1) x^{2(m-1)+1}/(2(m-1)+1) ! × (-1)x^2/[(2m)(2m+1)]
= a(m-1) ×　(-1)x^2/[(2m)(2m+1)]

これを利用したアルゴリズムです。

こんなかんじですね。

sx=0.0;
a=x;
n=1;
printf("sin(%lf) = %.8lf\n", xin, sin(x));
printf("----------------------------------------------\n");
do {
sx += a;
printf("%3d: sin(%lf) = %.8lf, a(%d)=%le\n", n, x, sx, n, a);
a = a*(-1)*x*x/(n+1)/(n+2);
n+=2;
} while(fabs(a)>1e-8);
printf("----------------------------------------------\n");
printf("a(%d)=%le\n", n, a);

xが大きな数になるとオーバーフローするので、
0<x<2πになるように前処理はしておいたほうがいいと思います。

No.3 回答者： tatsu99 回答日時： 2009/07/16 12:04

以下のようにして下さい。

-----------------------
#include <stdio.h>
// xのn乗を求める
double mypow(double x,int n)
{
int i;
double ans = x;
for (i = 1; i < n; i++){
ans *= x;
}
return ans;
}
// nの階乗を求める
double myfact(int n)
{
if (n == 1) return 1.0;
return (myfact(n-1) * (double)n);
}
int main(void)
{
int i;
int fugou = 1;
double x=-1;
double v,absv;
double sinx = 0.0;
printf("x=%f\n",x);
printf("sinx = %f\n",sinx);

for (i = 1; ; i += 2){
// xのi乗 / iの階乗 を求める
v = mypow(x,i) / myfact(i);
// 誤差の10の-8乗未満なら終了
//一旦、絶対値を求め、絶対値で判断する（その方がわかりやすい）
absv = v;
if (absv < 0.0) absv *= (-1.0); //負なら正にする（絶対値にする）
if (absv < 1.0e-8) break; //絶対値で比較
sinx += (v*fugou);
fugou *= (-1); //iが3,7,11...の時は引き算なのでマイナスにする
}
printf("sinx = %f\n",sinx);
return (0);
}
-----------------------

No.1 回答者： hitokotonusi 回答日時： 2009/07/16 10:08

＞a = (-1)*v*m/((n+1)*(n+2));

は

a = a*(-1)*x*x/(n+1)/(n+2);
(Cらしく書けば a *= (-1)*x*x/(n+1)/(n+2);)

です。n-2の時のaを掛け忘れています。変数v,mは不要ですが、
少しでも演算速度をあげるならm=x*xを【ループの外】に書いてx*xをmに置き換える。
nは整数のほうがいいでしょう。

>v=v*(-1);

で符号を変えた後、もう一度a=(-1)*・・・で-1を掛けているので、結局符号は変わっていませんね。

この回答へのお礼

間違えてた意味がわかりました。どうもありがとうございました。

お礼日時：2009/07/16 14:30