用 C 中的泰勒级数逼近 Sine(x) 并有很多问题

Question

我正在尝试使用泰勒级数和斯特林对阶乘的近似来近似 C 中的正弦 (x)，但是对于任何 n=>5，对于 n

泰勒级数

Stirling 的阶乘近似

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

int main(){

float x,n,s,i,e,p,f,r;
f=M_PI;
e=2.7182818;
s=0;

printf("What value of sine do you want to apporximate?");
    scanf("%f", &x);
printf("With what level of precision do you want to calculate it?");
    scanf("%f", &n);
for(i=0;i<=n; ++i);{
    r=((2*i)+1);
    p=(sqrt(2*r*f)*(pow((r/e),r)));
    s=s+(((pow((-1),i))*(pow(x,((2*i)+1))))/p);
}
printf("the value of sine at %f is %f",x,s);
}

Answer 1

这一行

    for(i = 0; i <= n; ++i);{

有一个额外的分号。你正在执行一个空循环。

Answer 2

根据您的公式，这是正确的代码，但它会生成错误的输出，因此您需要再次检查您的公式：

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

int main(){

double x,n,s,i,e,p,f;
f=M_PI;
e=2.7182818;
s=0;
int sign=0;// Adding this type to toggle the sign 

printf("What value of sine do you want to apporximate?");
    scanf("%lf", &x);// conversion specifier must be %lf for floating number
printf("With what level of precision do you want to calculate it?");
    scanf("%lf", &n);
for(i=1;i<=n; i=i+2){ // Correcting the for loop
    p=sqrt(2*i*f)*pow((i/e),i); 
    s=s+(pow(-1,sign++)*pow(x,i))/p;
}
printf("the value of sine at %f is %f",x,s);
}

Answer 3

这更容易兄弟

#include <iostream>
#include <cmath>
using namespace std;
double factorial(int X)
{
    double factorial = 1;
    for(int i=1; i<=X; i++)
    {
        factorial = factorial *i;
    }
    return factorial;
}
double Mysin(double x,double result)
{
        for(int i = 0;i<20;i++)
    {
        result+=pow((-1),i)*pow(x,((2*i)+1))/factorial((2*i)+1);
    }
    return result;
}
double Mycos(double x,double result)
{
    for(int i = 0;i<20;i++)
    {
        result+=pow(-1,i)*pow(x,2*i)/factorial(2*i);
    }
    return result;
}
double Mytan(double sine,double cosine)
{
    return sine/cosine;
}
double deg_to_rad(double x)
{
    double const pi = 3.14159265359;
    return x*pi/180;
}
int main()
{
    double x,result=0;
    cin>>x;
    cout<<"My sin: "<<Mysin(deg_to_rad(x),result)<<endl;
    cout<<"My cosin: "<<Mycos(deg_to_rad(x),result)<<endl;
    cout<<"My tan: "<<Mytan(Mysin(deg_to_rad(x),result),Mycos(deg_to_rad(x),result))<<endl;
    return 0;
}

Answer 4

由于sin() 是一个周期函数，我不应该超过一个周期来计算它。这简化了太多的数学运算，因为您永远不需要计算很大的阶乘数。实际上，您甚至不需要计算系列中每个项的阶乘，因为系数可以从最后一项导出，只需将前一个系数除以(n-1) 和n。如果您的输入仅限于一个周期（好吧，您不需要使用固定周期M_PI，您可以假设最大值为3.5，并通过仅减少除以M_PI 的模数。

一旦这样说，我们可以绑定你的最大误差，至于3.5 的最大输入，我们将3.5^n/n! 作为我们近似的最后一项，其中一些n 的边界小于一些修正我们需要计算的项数的最大误差。

我不会试图精确计算需要计算的项数，而是尝试做一些猜测，从推导算法和显示实际值（例如，3.2 的最大输入值）

这些是n 位置处术语的值，用于输入3.2 和

   n  | term at position n for input `3.2`
======+=================
   8  |  0.27269634
  12  |  0.00240693
  16  |  0.00000578
  18  |  0.00000019
  20  |  0.00000001
  21  |  7.9E-10

所以我们可以停止只计算系列的 20 项。 exp() 函数确实如此，它添加了所有术语并且是一个简单的函数。对于sin() 或cos()，如果您认为两者具有exp() 函数的相同项，您可以猜出更好的误差估计，（第一个只有奇数项，第二个只有偶数项)

 (x^n)/(n!) - (x^(n+2))/((n+2)!) = (n!*x^n*(1 - x^2/((n+1)*(n+2))))/n!

n > 3.2 表示每个词是

< x^n/n!

所以我们可以应用与指数相同的标准。

这意味着我们可以在某个点停下来......如果我们继续我们的表格，我们会看到，例如n > 30，总的累积期限小于5.3E-18，所以我们可以停在那里（例如double 号码，至少）。

#include <stdio.h>
#include <math.h> /* for the system sin() function */

double MySin(double x) /* x must be in the range [0..3.2] */
{
    int i;
    const int n = 30;
    double t = x, acum = x; /* first term, x/1! */
    x *= x; /* square the argument so we get x^2 in variable x */
    for (i = 3; i < n; i += 2) {
             t = -t * x / i / (i-1); /* mutiply by -1, x^2 and divide by i and (i-1) */
             acum += t; /* and add it to the accum */
    }
    return acum;
}

int main()
{
    double arg;
    for(;;) {
            if (scanf("%lg", &arg) != 1)
                    break;
            printf("MySin(%lg) = %lg; sin(%lg) = %lg\n",
                    arg, MySin(arg), arg, sin(arg));
    }
}

如果你利用 sin 函数的对称性，你可以将你的域减少到小于 1 的M_PI/4，你甚至可以在 18 次方处停止以获得大约 17 个有效数字（对于double) 让你更快犯罪。

最后，我们可以得到一个对所有真实域都有效的sin2()函数：

double sin2(double x)
{
    bool neg = false;
    int ip = x / 2.0 / M_PI;
    x -= 2.0 * M_PI * ip;                    /* reduce to first period [-2PI..2PI] */
    if (x < 0.0) x += 2.0*M_PI;              /* reduce to first period [0..2PI] */
    if (x > M_PI) { x -= M_PI; neg = true; } /* ... first period negative [ 0..PI ] */
    if (x > M_PI/2.0) x = M_PI - x;          /* reflection [0..PI/2] */
    return neg ? MySin(-x) : MySin(x);
}