时间戳

首页 > 解决方案 > 为什么我在 C 中的罪孽计算代码返回错误的值？

c - 为什么我在 C 中的罪孽计算代码返回错误的值？

问题描述

我想根据这个等式使用我自己的函数计算用户输入的正弦值：

sin(x) = sum_(i=0)^n (-1)^i * (x^(2 i + 1)/((2 i + 1)!))

方程的图片

我有这段代码，据我了解，我所做的与方程式中写的完全相同：

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

int faculty(int factor)
{
    int result = 1;

    if (factor > 0)
    {
        for (int i = 1; i <= factor; i++)
        {
            result = result * i;
        }
    }
    else
    {
        result = 1;
    }

    return result;
}

double my_sin(double x, int n)
{
    double my_sin = 0;
    for (int j = 0; j < n ; j++)
    {
        double i = (double)j;
        double faculty_j = (double)faculty(2*j+1);
        my_sin = my_sin + (pow((-1.0), i) * (pow(x, (double)(2.0 * i + 1.0)) / faculty_j));
    }

    return my_sin;
}

int main()
{
    int n = 0;
    double x = 0;

    printf("x=");
    scanf("%lf", &x);

    printf("n=");
    scanf("%i", &n);

    printf("sin(%i)=", (int)x);
    printf("%lf\n", my_sin(x, n));

    return 0;
}

但是，例如，当我使用时x = 8，n = 5我得到sin(8)=149了结果。我现在尝试调试代码一段时间，但我不知道问题可能出在哪里或如何找出问题所在。

更新代码：

#include <stdio.h>

long int factorial(int factor)
{
    long int result = 1;

    if (factor > 0)
    {
        for (int i = 1; i <= factor; i++)
        {
            result = result * i;
        }
    }
    else
    {
        result = 1;
    }

    return result;
}

double my_pow(double a, double b)
{
    if (b == 0)
    {
        return 1;
    }

    double result = a;
    double increment = a;
    double i, j;

    for (i = 1; i < b; i++)
    {
        for (j = 1; j < a; j++)
        {
            result += increment;
        }
        increment = result;
    }
    return result;
}

double my_sin(double x, int n)
{
    double my_sin = 0;
    for (int j = 0; j < n ; j++)
    {
        double i = (double)j;
        double faculty_j = (double)factorial(2*i+1);
        my_sin = my_sin + (my_pow((-1.0), i) * (my_pow(x, 2.0 * i + 1.0) / faculty_j));
    }

    return my_sin;
}

int main()
{
    int n = 0;
    double x = 0;

    printf("x=");
    scanf_s("%lf", &x);

    printf("n=");
    scanf_s("%i", &n);

    printf("sin(%i)=", (int)x);
    printf("%lf\n", my_sin(x, n));

    return 0;
}

标签： cequation-solving

解决方案

这里的问题不一定是你的代码。这只是您对方程式的期望。我在 Wolfram Alpha 中测试了 n=5 的泰勒级数，并使用查询从级数中减去了 sin(x)

(sum (-1)^i * x^(2i+1)/(2i+1)!, i=0 to 5) - sin(x)

查看错误。

https://www.wolframalpha.com/input/?i=%28sum+%28-1%29%5Ei+*+x%5E%282i%2B1%29%2F%282i%2B1%29%21%2C+i %3D0+到+5%29+-+sin%28x%29

在此处输入图像描述

如您所见，当 x 大约为 0 时，该级数的近似值非常好，但当 x 超过 3 时，误差开始变得很大，并且增长得非常快。在 x=8 时，错误的大小为 67，这显然是无用的。要获得合理的近似值，您需要使用 n=15 左右。只有五个太少了。

您应该在这里做的是利用 sin(x) = sin(x + k * 2 * PI) 的事实，其中 k 是任意整数。如果 x 是负数，你可以简单地使用 sin(x) = -sin(-x) 这个事实。我正在展示一个简单的例子来说明它是如何完成的：

double my_sin(double x, int n)
{
    double my_sin = 0;
    double sign = 1.0;

    // sin(-x) = -sin(x)
    // This does not improve accuracy. It's just to make the rest simpler
    if(x<0) {
        x=-x;
        sign *= -1;
    }

    // sin(x) = sin(k*2*PI + x)
    x -= 2*PI*floor(x/(2*PI));

    // Continue as usual

    return sign * my_sin;
}

您希望尽可能接近 x=0。为了进一步改进这一点，您可以使用 sin(x) = -sin(x+PI) 的事实，如下所示。

以上足以将所有内容保持在区间 [0, PI] 中，但实际上我们可以通过使用一些更聪明的数学来做得更好。我们可以使用 sin(x) 关于 x=PI/2 + k*PI 对称的事实。所以 sin(PI/2 -x) = sin(PI/2 + x)。利用它将使您保持在区间 [0, PI/2] 中，当您处于该区间时，n=2 足以使误差低于 0.005。对于 n=5，误差低于 10^-7。您可以在前面的步骤之后添加它：

// sin(x) = -sin(x+PI)
if(x>PI) {
    x-=PI;
    sign *= -1;
}

// sin(PI/2 -x) = sin(PI/2 + x)
if(x>PI/2) 
    x = PI - x;

但是，您的代码中有一个错误。请记住，应该包括案例 j=n，因此请更改

for (int j = 0; j < n ; j++)

至

for (int j = 0; j <= n ; j++)

另一种方法是始终在参数中添加 1 来调用函数。

只是为了清楚起见。请记住，泰勒级数期望x以弧度为单位。如果你使用学位，你会得到错误的结果。

为了完整起见，这里是完整的功能和一些其他不影响正确性但性能和可读性的修复：

double my_sin(double x, int n)
{
    double ret = 0;
    double sign = 1.0;

    if(x<0) {
        x=-x;
        sign *= -1;
    }

    x -= 2*PI*floor(x/(2*PI));

    if(x>PI) {
        x-=PI;
        sign *= -1;
    }

    if(x>PI/2) 
        x = PI - x;

    size_t denominator = 1;
    double numerator = x;
    int s = -1;

    for (int j = 0; j <= n ; j++) {
        denominator *= 2*j + 1;
        s *= -1;

        ret += s * numerator / denominator;

        numerator *= x*x;
    }

    return sign*ret;
}

以一定的精度为目标

假设您想要 3 个正确的十进制数字。然后我们可以利用分子比分母增长慢的事实，做这样的事情：

    size_t denominator = 1;
    double numerator = x;
    double term;
    int s = -1;
    int j = 0;
    double tolerance = 0.001;

    do {
        denominator *= 2*j + 1;
        j++;
        s *= -1;

        term = numerator / denominator;

        ret += s * term;

        numerator *= x*x;
    } while(abs(term) > tolerance);

推荐阅读