正弦和余数的近似值

Approximating of sine and the remainder

本文关键字：近似值 余数      更新时间：2023-10-16

练习的目标是使用麦克劳林级数公式评估sin。

#include <stdio.h>
#include <math.h>
double factorial(int n);
int main(void) {
    double x, p, r;
    int i, n;
    printf("Enter a positive double number and a non - negative integer : n");
    scanf("%lf%d", &x, &n);
    if (x <= 0) {
        printf("Error: first argument must be a positive double.n");
        return -1;
    }
    if (n < 0) {
        printf("Error: second argument must be a non - negative integer.n");
        return -1;
    }
    p = 0;
    for (i = 0; i <= n; i++)
    {
        p += pow(-1, i) / factorial(2 * i + 1) * pow(x, 2 * i + 1);
    }
    r = fabs(sin(x) - p);
    printf("The %d-th order Maclaurin polynomial function at x=%f is %f, with an error approximation of %f.n", n, x, p, r);

    getch();
    return 0;
}
double factorial(int n)
{
    int i;
    long result = 1;
    for (i = 1; i <= n; i++)
        result *= i;
    return result;
}

我得到输入"12 16"的奇怪结果。为什么？

这里有三个问题。

  • 正如Mukit Chowdhury回答的那样，longs不能容纳大的阶乘。但是，正如GregS评论的那样，16！应该没有问题。当你使用多头来保持 21 时，你应该期待奇怪的结果！或更大，但您无需更改此输入的阶乘函数。不过，您可能应该改用类似以下内容的内容：

    double factorial(int n)
{
    int i;
    double result = 1;
    for (i = 1; i <= n; i++)
        result *= i;
    return result;
}

  • 在输入"12 16"上，您的代码声称正在计算 16 阶麦克劳林多项式，但它计算的是 33 阶麦克劳林多项式。16阶多项式的项数高达 -x^15/15！+ 0x^16。解决此问题的一种方法是更正 for 循环，如下所示：

    for (i = 1; i <= n; i+=2)
{
    p += pow(-1, (i-1)/2) / factorial(i) * pow(x, i);
}

    因此，您的代码遇到了阶乘问题，但这只是因为您正在计算额外的项。如果计算的项最多为 -x^15/15！，您应该能够正确计算多项式的值。

    • 16
  • 阶麦克劳林多项式在12处的实际值是-4306.756...这可能不是你所期望的，这可能是练习的一部分。为了获得准确的近似值，您应该期望最后一个项很小，因此您需要 n！超过 x^n。根据斯特林的近似，n！~ (n/e(^n，所以你需要 n> e*x，其中 e = 2.71828...，所以 n>=33。此时，误差为 0.005，将 n 增加 c 会使误差的大小减小大约 e^c 的因子。
  • 当您在产生小的最终结果的过程中减去大数字时，您应该预料到双精度算术中的大误差。这可能不是问题，因为最大的项只有大约 2^14 的量级。你仍然会得到足够的精度，以至于你不会注意到你不能通过添加更多的术语来接近sin(12(。

麦克劳林或泰勒级数的两个项之间的商为

-x*x/(2*k*(2k+1))

可以有利地利用这一点来避免所有幂和阶乘及其溢出。

mxx = -x*x;
term = mxx / 6;
sum = 1+term;
k=2;
while( not good enough )
    term = term*mxx/(2*k*(2*k+1));
    sum = sum + term;
    k = k+1;
return sum*x
在这里我会

做很多不同的事情。

我

有点偏向于编写数字软件，该软件不仅可以获得正确的结果，而且可以快速获得结果，但是我在这里看到了很多浪费的计算。例如，考虑序列的两个连续非零项(x^13(/(13！( 和 -(x^15(/(15！(。如果您已经知道(x^13(/(13！(，你需要做多少计算才能得到 -(x^15(/(15！(？答案是，比计算(x^13(/(13！(所花费的时间要少得多。首先。如果你奴性地遵循通常的公式对于麦克劳林级数并重新计算每个新项的阶乘，为了得到15！您将重复已经完成的所有计算对于 13！，然后只执行两个新的乘法。在不浪费计算和没有浪费计算的情况下有效地计算此序列引入不必要的大数字(以及所有可能的问题即使您对所有内容都使用浮点数，它们也可能导致(，看看一个非零项之间的比率和下一个。它易于计算，不涉及拖入pow函数。一个好的算法会让你轻松增加术语的数量，直到最后一个术语是浮点数的精度顺序，前提是您使用合理的输入值 x。(x不需要大于2 * pi，或者真的没有充分的理由将 x 设置为大于 pi/2，因为任何更大的正弦x 的值可以从该范围内输入的正弦值中找到。

以下是我在 Java 中的操作(易于转换为 C++(：

package math;
/**
 * Sine using Maclauren series
 * @author Michael
 * @link http://stackoverflow.com/questions/29445615/approximating-of-sine-and-the-remainder
 * @since 4/4/2015 8:37 PM
 */
public class Sine {
    private static final double TWO_PI = 2.0*Math.PI;
    private static final int numPoints = 21;
    private static final int numTerms = 21;
    public static void main(String[] args) {
        double x = -Math.PI;
        double dx = 2.0*Math.PI/(numPoints-1);
        for (int i = 0; i < numPoints; ++i) {
            System.out.println(String.format("# terms: %d angle (radians) %10.6f sine: %15.10f", numTerms, x, sine(x, numTerms)));
            x += dx;
        }
    }
    public static double sine(double radians, int numTerms) {
        double value = 0.0;
        // Start by making sure the angle -pi/2 <= x <= +pi/2
        double x = normalizeAngle(radians);
        double term = x;
        for (int n = 3; n < numTerms; n += 2) {
            value += term;
            term *= -x*x/n/(n-1);
        }
        return value;
    }
    public static double normalizeAngle(double radians) {
        return radians - TWO_PI*Math.floor((radians+Math.PI)/TWO_PI);
    }
}

下面是输出：

java math.Sine
# terms: 21 angle (radians)  -3.141593 sine:   -0.0000000224
# terms: 21 angle (radians)  -2.827433 sine:   -0.3090169974
# terms: 21 angle (radians)  -2.513274 sine:   -0.5877852526
# terms: 21 angle (radians)  -2.199115 sine:   -0.8090169944
# terms: 21 angle (radians)  -1.884956 sine:   -0.9510565163
# terms: 21 angle (radians)  -1.570796 sine:   -1.0000000000
# terms: 21 angle (radians)  -1.256637 sine:   -0.9510565163
# terms: 21 angle (radians)  -0.942478 sine:   -0.8090169944
# terms: 21 angle (radians)  -0.628319 sine:   -0.5877852523
# terms: 21 angle (radians)  -0.314159 sine:   -0.3090169944
# terms: 21 angle (radians)   0.000000 sine:    0.0000000000
# terms: 21 angle (radians)   0.314159 sine:    0.3090169944
# terms: 21 angle (radians)   0.628319 sine:    0.5877852523
# terms: 21 angle (radians)   0.942478 sine:    0.8090169944
# terms: 21 angle (radians)   1.256637 sine:    0.9510565163
# terms: 21 angle (radians)   1.570796 sine:    1.0000000000
# terms: 21 angle (radians)   1.884956 sine:    0.9510565163
# terms: 21 angle (radians)   2.199115 sine:    0.8090169944
# terms: 21 angle (radians)   2.513274 sine:    0.5877852526
# terms: 21 angle (radians)   2.827433 sine:    0.3090169974
# terms: 21 angle (radians)   3.141593 sine:   -0.0000000224
Process finished with exit code 0

你的阶乘((返回一个长，而它的返回类型是双精度。

最重要的是，长不能包含16！哪个双可以。

相关文章：