【问题标题】:taylor series with error at most 10^-3泰勒级数最多有 10^-3 的误差
【发布时间】:2015-12-10 02:54:46
【问题描述】:

我正在尝试计算cos(x) 的泰勒级数,最多有错误10^-3 和所有x ∈ [-pi/4, pi/4],这意味着我的错误需要小于0.001。我可以修改 for 循环中的 x += 以获得不同的结果。我尝试了几个数字,但它从来没有变成小于 0.001 的错误。

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

float cosine(float x, int j)
{
    float val = 1;
    for (int k = j - 1; k >= 0; --k)
        val = 1 - x*x/(2*k+2)/(2*k+1)*val;
    return val;
}

int main( void )
{
   for( double x = 0; x <= PI/4; x += 0.9999 )
   {

       if(cosine(x, 2) <= 0.001)
       {
           printf("cos(x) : %10g    %10g    %10g\n", x, cos(x), cosine(x, 2));
       }
       printf("cos(x) : %10g    %10g    %10g\n", x, cos(x), cosine(x, 2));
    }

    return 0;
}

我也在为e^x 做这件事。对于这部分,x must in [-2,2]

float exponential(int n, float x)
{
    float sum = 1.0f; // initialize sum of series

    for (int i = n - 1; i > 0; --i )
        sum = 1 + x * sum / i;

    return sum;
}

int main( void )
{
    // change the number of x in for loop so you can have different range
    for( float x = -2.0f; x <= 2.0f; x += 1.587 )
    {
        // change the frist parameter to have different n value
        if(exponential(5, x) <= 0.001)
        {
            printf("e^x = %f\n", exponential(5, x));
        }
    printf("e^x = %f\n", exponential(5, x));
    }

    return 0;
}

但是每当我更改 for 循环中的项数时,它总是会出现大于 1 的错误。我应该如何将其更改为小于 10^-3 的错误?

谢谢!

【问题讨论】:

  • 如果将exponential 循环中的计算更改为sum = 1.0 + x * sum / (float) i 会发生什么?
  • 好的,我看到它应该是浮动的,但除此之外,没有什么真正改变,输出保持不变。
  • 在您检查值if(cosine(x, 2) &lt;= 0.001)if(exponential(5, x) &lt;= 0.001) 的这些行中,您不应该与cosexp 的TRUE 值进行比较吗?您在这里检查的似乎是函数是否接近于零。我觉得应该是if (abs(cosine(x,2) - cos(x)) &gt; .0001) { ... }

标签: c taylor-series function-approximation


【解决方案1】:

我的理解是,为了提高精度,您需要考虑泰勒级数中的更多项。例如,考虑一下当 您尝试通过泰勒级数计算 e(1)。

$e(x) = \sum\limits_{n=0}^{\infty} frac{x^n}{n!}$

我们可以考虑e(1)展开式中的前几项:

n             value of nth term           sum
0        x^0/0! = 1                   1
1        x^1/1! = 1                   2
2        x^2/2! = 0.5                 2.5
3        x^3/3! = 0.16667             2.66667
4        x^4/4! = 0.04167             2.70834

您应该注意到两件事,首先,随着我们添加更多项,我们越来越接近 e(1) 的确切值,而且连续和之间的差异越来越小。

所以,e(x) 的实现可以写成:

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

typedef float (*term)(int, int);
float evalSum(int, int, int, term);
float expTerm(int, int);
int fact(int);
int mypow(int, int);
bool sgn(float);

const int maxTerm = 10;         // number of terms to evaluate in series
const float epsilon = 0.001;    // the accepted error

int main(void)
{
    // change these values to modify the range and increment 
    float   start = -2;
    float   end = 2;
    float   inc = 1;

    for(int x = start; x <= end; x += inc)
    {
        float value = 0;
        float prev = 0;

        for(int ndx = 0; ndx < maxTerm; ndx++)
        {
            value = evalSum(0, ndx, x, expTerm);

            float diff = fabs(value-prev);
            if((sgn(value) && sgn(prev)) && (diff < epsilon))
                 break;
            else
                 prev = value;
        }

        printf("the approximate value of exp(%d) is %f\n", x, value);
    }

    return 0;
}

我曾经猜测我们不需要在扩展中使用超过十个项来达到所需的精度,因此内部 for 循环是我们在范围内循环 n 的值0,10]。

此外,我们有几行专门用于检查我们是否达到了所需的精度。首先我计算当前评价与上一次评价的差值的绝对值,取绝对差值。检查差值是否小于我们的 epsilon 值 (1E-3) 是提前退出循环的标准之一。我还需要检查当前值和之前值的符号是否相同,因为计算 e(-1) 的值时存在一些波动,这就是条件中的第一个子句正在做的事情。

float evalSum(int start, int end, int val, term fnct)
{
    float sum = 0;
    for(int n = start; n <= end; n++)
    {
        sum += fnct(n, val);
    }

   return sum;
}

这是我编写的一个效用函数,用于评估系列的前 n 项。 start 是起始值(此代码始终为 0),end 是结束值。最后一个参数是一个指向函数的指针,该函数表示如何计算给定项。在此代码中,fnct 可以是指向任何接受整数参数并返回浮点数的函数的指针。

float expTerm(int n, int x)
{
    return (float)mypow(x,n)/(float)fact(n);
}

隐藏在这个单行函数中的是大部分工作发生的地方。此函数表示 e(n) 的泰勒展开的封闭形式。仔细查看上面的内容,您应该能够看到我们正在为给定的 x 和 n 值计算 $\fract{x^n}{n!}$。作为提示,为了执行余弦部分,您需要创建一个函数来评估 cos 的泰勒展开中的一项的闭项。这是由 $(-1)^n\fact{x^{2n}}{(2n)!}$ 给出的。

int fact(int n)
{
    if(0 == n)
        return 1;             // by defination
    else if(1 == n)
        return 1;
    else
        return n*fact(n-1);
}

这只是阶乘函数的标准实现。这里没什么特别的。

int mypow(int base, int exp)
{
    int result = 1;

    while(exp)
    {
        if(exp&1)              // b&1 quick check for odd power
        {
            result *= base;
        }

        exp >>=1;              // exp >>= 1 quick division by 2
        base *= base;
    }

    return result;
}

一个用于求幂的自定义函数。我们当然可以使用来自&lt;math.h&gt; 的版本,但是因为我知道我们只会做整数幂,所以我们可以编写一个优化的版本。提示:在计算余弦时,您可能需要使用 &lt;math.h&gt; 中的版本来处理浮点基数。

bool sgn(float x)
{
    if(x < 0) return false;
    else return true;
}

一个非常简单的函数,用于确定浮点值的符号,返回真为正,否则为假。

这段代码是在我的 Ubuntu-14.04 上使用 gcc 版本 4.8.4 编译的:

******@crossbow:~/personal/projects$ gcc -std=c99 -pedantic -Wall series.c -o series
******@crossbow:~/personal/projects$ ./series
the approximate value of exp(-2) is 0.135097
the approximate value of exp(-1) is 0.367857
the approximate value of exp(0) is 1.000000
the approximate value of exp(1) is 2.718254
the approximate value of exp(2) is 7.388713

使用 bc 给出的预期值是:

******@crossbow:~$ bc -l
bc 1.06.95
Copyright 1991-1994, 1997, 1998, 2000, 2004, 2006 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
e(-2)
.13533528323661269189
e(-1)
.36787944117144232159
e(0)
1.00000000000000000000
e(1)
2.71828182845904523536
e(2)
7.38905609893065022723

如您所见,这些值完全在您要求的公差范围内。我把它作为一个练习来做余弦部分。

希望这会有所帮助,
-T

【讨论】:

    【解决方案2】:

    expcos 具有在实线上处处收敛的幂级数。对于任何有界区间,例如[-pi/4, pi/4][-2, 2],幂级数不仅逐点收敛,而且一致地收敛到 expcos

    逐点收敛意味着对于该区域中的任何x,以及任何epsilon &gt; 0,您都可以选择一个足够大的N,以便您从泰勒级数的第一个N 项中得到的近似值在epsilon 的真实值。但是,通过逐点收敛,N 对于某些 x 可能很小,而对于其他 x 可能很大,并且由于 x 的数量是无限的,因此可能没有有限的 N 可以容纳它们。对于某些功能,有时确实会发生这种情况。

    均匀收敛意味着对于任何epsilon &gt; 0,您可以选择一个足够大的N,以便该区域中每个x 的近似值都在epsilon 之内。这就是您正在寻找的近似值,并且可以保证这就是您所拥有的那种收敛。

    原则上,您可以查看expcos 在任何有限域上一致收敛的证明之一,坐下来说“如果我们采用epsilon = .001 会怎样,并且区域要... ",并使用笔和纸计算N 上的一些有限界限。然而,这些证明中的大多数将在某些步骤中使用一些不明确的估计,因此您计算的 N 的值将大于必要的值——可能要大得多。将N 作为一个变量来实现它会更简单,然后像在代码中那样使用for循环检查值,看看你必须使它有多大才能使错误小于@987654346 @无处不在。

    所以,我无法确定您需要选择的 N 的正确值是多少,但数学可以保证,如果您继续尝试更大的值,最终您会找到一个有效的值。

    【讨论】:

      猜你喜欢
      • 2015-11-23
      • 2021-12-26
      • 2014-04-03
      • 1970-01-01
      • 2017-12-02
      • 2014-02-23
      • 1970-01-01
      • 1970-01-01
      • 2011-12-06
      相关资源
      最近更新 更多