5

对于我的研究，我必须编写一个算法来sin()使用这个函数进行计算：

但是，在我的算法中，我必须将 X 的值保持在 0 和 Pi/2 之间。所以，我写了我的算法，但所有的结果都是错误的。

这是我的代码：

double sinX(double x){
    double resultat = 0;
    int i;
    if(x < 0 || x > M_PI_2)
        x = fmod(x,M_PI_2);

    for(i = 1;i<=30;i++){
        resultat += -1 * ((x*x)/(2*i*(2*i+1)))*(pow(-1,i-1))*((pow(x,2*i-1))/(factorielle(2*i-1)));
    }

    return resultat;
}

我没找到原因。你能帮助我吗？

这里是 X 的几个值和 fmod 的结果

1 / 1
2 / 0.429204
3 / 1.4292
4 / 0.858407
5 / 0.287611
6 / 1.28761
7 / 0.716815
8 / 0.146018
9 / 1.14602
10 / 0.575222
11 / 0.00442571
12 / 1.00443
13 / 0.433629
14 / 1.43363
15 / 0.862833
16 / 0.292037
17 / 1.29204
18 / 0.72124
19 / 0.150444
20 / 1.15044

和算法的结果

1 / -0.158529
2 / -0.0130568
3 / -0.439211
4 / -0.101605
5 / -0.00394883
6 / -0.327441
7 / -0.0598281
8 / -0.000518332
9 / -0.234888
10 / -0.0312009
11 / -1.44477e-008
12 / -0.160572
13 / -0.0134623
14 / -0.443022
15 / -0.103145
16 / -0.00413342
17 / -0.330639
18 / -0.0609237
19 / -0.000566869
20 / -0.237499

这是我的“工厂”定义

double factorielle(double x){
    double resultat = 1;
    int i;

    if(x != 0){
        for (i=2;i<=x;i++)
        {
            resultat *= i;
        }
    }
    else{
        resultat = 1;
    }

    return resultat;
}

和价值观：

1 / 1
2 / 2
3 / 6
4 / 24
5 / 120
6 / 720
7 / 5040
8 / 40320
9 / 362880
10 / 3.6288e+006
11 / 3.99168e+007
12 / 4.79002e+008
13 / 6.22702e+009
14 / 8.71783e+010
15 / 1.30767e+012
16 / 2.09228e+013
17 / 3.55687e+014
18 / 6.40237e+015
19 / 1.21645e+017
20 / 2.4329e+018
4

4 回答 4

8

你误解了你展示的第二个公式的目的。这个想法是您使用该公式来计算前一项的总和中的每个项，从而使您无需使用任何pow或factorial调用。

#include <stdio.h>

double sinX(double x) {
  double term, total_so_far;
  int i;

  term = x;  /* First term in the expansion. */
  total_so_far = 0.0;
  for (i = 1; i <= 30; i++) {
    /* Add current term to sum. */
    total_so_far += term;
    /* Compute next term from the current one. */
    term *= -(x * x) / (2*i) / (2*i + 1);
  }
  return total_so_far;
}


int main(void) {
  /* testing */
  double x;
  int i;

  for (i = 0; i <= 10; i++) {
    x = i / 10.0;
    printf("sin(%f) is %f\n", x, sinX(x));
  }
  return 0;
}

以及在我的机器上运行此代码的结果：

sin(0.000000) is 0.000000
sin(0.100000) is 0.099833
sin(0.200000) is 0.198669
sin(0.300000) is 0.295520
sin(0.400000) is 0.389418
sin(0.500000) is 0.479426
sin(0.600000) is 0.564642
sin(0.700000) is 0.644218
sin(0.800000) is 0.717356
sin(0.900000) is 0.783327
sin(1.000000) is 0.841471

这应该给你合理的结果范围0为pi / 2。在该范围之外，您需要对正在使用的减少更加聪明：简单地减少模pi / 2不会给出正确的结果。（提示：减少 modulo 是安全的2 * pi，因为sin函数是周期性的 period 2 * pi。现在使用sin函数的对称性来减少到 的范围0。pi / 2）

编辑对当前代码给出错误结果的原因的解释：除了有缺陷的归约步骤之外，在您的总和中，您从 term 开始i = 1。但是第一个词应该是for i = 0（就是x词，i=1词就是-x^3 / 3!词）。一个快速而肮脏的解决方法是删除归约步骤，并将resultat变量初始化为x而不是0. 那应该会给您带来好的结果 small x，然后您就可以弄清楚如何替换减少步骤。不过，如果您真的打算使用显式阶乘和幂调用来计算答案，我会感到惊讶 - 我几乎可以肯定，您应该如上所述从前一项计算每个项。

于 2013-11-11T16:01:21.977 回答
3

您的代码有两个问题：

  1. sin(x+k*π/2)不一定等于sin(x)
  2. 你对这个词的表达有点混乱。说明似乎建议您从上一个术语计算系列中的下一个术语。从值开始，i=0然后使用问题中的等式计算每次迭代中的下一项。
于 2013-11-11T15:59:40.853 回答
0

最后，我听从了你的指示。这是我的最终代码：

double sinX(double x)
{
    double result = 1.0;
    double term_i = 1.0;
    int i = 2;

    x = fmod(x, 2*M_PI);

    for(i = 2; i<= 30; i+=2)
    {
        term_i = (-term_i * (x*x)) / (i*(i+1));
        result += term_i;
    }

    return x * result;
}
于 2013-11-11T19:39:25.657 回答
0

关于 OP 发布答案的术语数量的想法。

只要先执行一些范围限制，例如fmod()，就可以合理地动态确定所需的项数。（对 x 使用 1 到 23 次迭代：0 到 2*pi。）

double sinX1(double x)
{
    double result = 1.0;
    double term_i = 1.0;
    int i = 2;

    x = fmod(x, 2*M_PI);

    // for(i = 2; i<= 30; i+=2)
    for(i = 2; ((1.0 + term_i) != 1.0); i+=2)
    {
        term_i = (-term_i * (x*x)) / (i*(i+1));
        result += term_i;
    }
    return x * result;
}
于 2013-11-13T16:05:37.540 回答