（六）算法 -- 4. 级数展开

4. 级数展开

• 级数展开(series expansion)
  函数的值通过对级数项的求和来估计。如果每加入一个新的项都会使总和更接近于一个所期望的值，则级数是收敛的，可以用此级数估计结果。

4.1 Zeno悖论

• Zeno悖论(Zeno’s paradox)
  横穿房间这个过程中的每一步都是走到剩余距离的一半。

已经走过的每段距离的总和，可以表示为下式：
$$\frac{1}{2}+\frac{1}{4}+\frac{1}{8}+\frac{1}{16}+\frac{1}{32}+\frac{1}{64}+\frac{1}{128}+...$$

如果从1开始给这些项编号，则和式中的第$i$项可写为：

$$\frac{1}{2^i}$$

第$i$项涉及$i$次幂的数学级数称为幂级数（power series）。

Zeno级数有无限项，但有有限和（为1）。

计算Zeno级数和，代码实现如下：

#include <stdio.h>

/*  Calcualte the sum of Zeno's series. */

/* Main Program */
main() {
    double sum, term;
    sum = 0.0;
    term = 0.5;

    while (sum != sum + term) {
        sum += term;
        term /= 2;
    }
    printf("The sum of Zeno's series is: %g\n", sum);
}

4.2 泰勒级数近似的实现

需要写一个函数，通过幂级数展开计算平方根的近似值：

$$\sqrt{x}\cong 1+\frac{1}{2}(x-1)-\frac{1}{4}\frac{(x-1{)}^2}{2!}+\frac{3}{8}\frac{(x-1{)}^3}{3!}-\frac{15}{16}\frac{(x-1{)}^4}{4!}+...$$

最有效的策略是找出用前一项计算新的项的方法。

每一项都有三个独立的部分：
$$corff\frac{xpower}{factorial}$$

对于xpower，将当前的xpower乘以x-1就得到了下一项的xpower。

对于factorial，第i项的值总是i!。任何数i的阶乘是i乘以比它小的那个数的阶乘。
为了准备第(i+1)个循环周期，需要把当前的factorial乘以i+1。

到第(i+1)个循环周期的coeff值，需要把第i个周期的值乘以$\frac{1}{2}$ - i。

代码实现如下：

#include <stdio.h>

/* Estimate square root using Taylor expansion when 0<x<2. */


/* Function Prototype */
double TSqrt(double x);


/* Main Program */
main() {
    double x;
    printf("Input a number x (0<x<2): ");
    scanf("%lf", &x);
    if (x <= 0 | x >= 2) {
        printf("Number x should meet the condition: 0<x<2.\n");
        exit(0);
    }
    printf("Estimated square root is: %g\n", TSqrt(x));
}


/* Function */
double TSqrt(double x)
{
    double corff, xpower, factorial, sum, term;
    int i;
    sum = 0;
    corff = xpower = factorial = term = 1;

    for (i = 0; sum != sum + term; i++) {
        sum += term;
        corff *= (0.5 - i);
        xpower *= (x - 1);
        factorial *= (i + 1);
        term = corff * xpower / factorial;
    }
    return sum;
}

TSqrt函数作为泰勒级数展开，平方根函数的公式仅当参数值位于一个有限范围内时才有效，在该范围内计算趋于收敛。该范围被称为收敛半径(radius of convergence)。

泰勒级数公式希望x处于范围0<x<2。

参考
《C语言的科学和艺术》 —— 第6章 算法