0

I have been working on a program for Taylor Series and used long double as the number format to allow calculation of large numbers. My program works pretty fine for positive exponents, yet it fails when it comes to negative exponents. The problem is that I am getting very large positive number when I calculate exp(-x) for some x. What may be the reason behind this? Thanks for your help in advance. You can see my code here: 

#include <stdio.h>
#include <math.h>
//We need to write a factorial function beforehand, since we
//have factorial in the denominators.
//Remembering that factorials are defined for integers; it is
//possible to define factorials of non-integer numbers using
//Gamma Function but we will omit that.
//We first declare the factorial function as follows:
long double factorial (double);
//Long long integer format only allows numbers in the order of 10^18 so 
//we shall use the sign bit in order to increase our range.
//Now we define it,
long double
factorial(double n)
{
//Here s is the free parameter which is increased by one in each step and
//pro is the initial product and by setting pro to be 0 we also cover the
//case of zero factorial.
    int s = 1;
    long double pro = 1;
    if (n < 0)
        printf("Factorial is not defined for a negative number \n");
    else {
    while (n >= s) { 
    pro *= s;
    s++;
    }
    return pro;
    }
}

int main ()
{
    long double x[13] = { 1, 5, 10, 15, 20, 50, 100, -1, -5, -10, -20, -50, -100};
//Here an array named "calc" is defined to store 
//the values of x.

//int k;
////The upper index controls the accuracy of the Taylor Series, so
////it is suitable to make it an adjustable parameter. 
int p = 135;
long double series[13][p];
long double sum = 0;
int i, k;
for (i = 0; i <= 12;i++) {
for (k = 0; k <= p; k++){
    series[i][k] = pow(x[i], k)/( factorial(k));
    sum += series[i][k];
}
printf("Approximation for x = %Lf is %Lf \n", x[i], sum);
}
printf("%Lf \n", factorial(100));
}
4

3 回答 3

5

这只是数值分析的数学主题。MacLaurin 级数e^x为所有人收敛x，但让我们看看为什么它对e^(-10).

e^x = 1 + x + x^2/2 + x^3/6 + x^4/24 + x^5/120 + x^6/720 + x^7/5040 + ... +x^n/n! + ...

e^(-10) = 1 - 10 + 100/2 - 1000/6 + 10000/24 -100000/120 + ...

该系列中最大的术语是什么？10^10/10!, 大约是2755.7319224. e^(-10) 大约的真实值是多少0.00004539992。将系列相加会丢失 9 位精度，而您根本没有。

如果你找到e^(10)并采取了倒数，你会相当安全。如果您通过将 (1/e) 乘以 10 次直接计算 e^(-10)，那么您也是安全的。但是，任何具有比真实答案大得多的交替项的序列都会导致这些问题。

即使对于范围有限的函数，MacLaurin 级数也不会在实践中使用。例如，首先采用 trig 函数的参数，并使用周期性和 trig 恒等式将参数减少到区间0 < θ < π/4。然后，人们经常应用切比雪夫近似来均匀地减少误差。在其他情况下，连分数和 Pade 近似值优于三角级数。贝塞尔函数最好通过向后递归来完成。

看一本好的数值分析书。Forman Acton 的“通常有效的数值方法”是老式的，但很好。

于 2013-10-07T22:03:00.040 回答
2

你没有把你的总和归零。您只是将每个新结果添加到前一个结果中。

添加为第一个循环sum = 0;中的第一个语句。for

道德：总是为应该是函数的函数创建函数。在这种情况下，编写一个单独的 exp_taylor() 函数，而不仅仅是内联编写它。

于 2013-10-07T21:50:23.387 回答
0

在计算每个系列后，您需要将总和重置为零：

int main ()
{
  long double x[13] = { 1, 5, 10, 15, 20, 50, 100, -1, -5, -10, -20, -50, -100}; 
  const int p = 135;
  long double series[13][p];
  int i, k;
  for (i = 0; i <= 12;i++) {
    long double sum = 0;  // LOOK HERE
    for (k = 0; k <= p; k++){
      series[i][k] = pow(x[i], k)/( factorial(k));
      sum += series[i][k];
    }
    printf("Approximation for x = %Lf is %Lf \n", x[i], sum);
  }
}

另请注意，您正在扩展大约 x=0 并且可以预期 x 远离 0 的显着错误。

于 2013-10-07T22:10:03.570 回答