用泰勒级数估算平方根需展开的阶数

要估算一个数的平方根，可以用泰勒级数，而需要展开的阶数和所需的精度以及求根的数有关。

多彩电脑

1700人浏览 · 2024-09-17 17:44:30

多彩电脑 · 2024-09-17 17:44:30 发布

我们知道，一个可导函数，可以用泰勒级数展开，这是一种估算函数值的方法。泰勒级数展开的基本公式是[1]：

$f(x)={f(x_0)\over 0!}+{{df(x_0)\over dx}\over 1!}(x-x_0)+{{d^2f(x_0)\over dx^2}\over 2!}(x-x_0)^2+. . . +{{d^nf(x_0)\over dx^n}\over n!}(x-x_0)^n+. . .$

用泰勒级数估算函数值的方法是只取前几项，忽略后面的无限多项。但泰勒级数并非必然收敛[2]，所以需要凑出能使泰勒级数收敛的 $x_0$ ，且 ${d^nf(x_0)\over dx^n}$ 的计算要容易，才能准确地用泰勒级数估算函数值。

博文[3]介绍了一种用泰勒级数估算平方根的方法，大致思路就是将平方根转化为 $n\sqrt{x}$ ，其中 $n$ 为正整数， $x$ 是一个接近1的数。但该文的方法只取泰勒展开式的二阶近似，有时准确，有时不准确。当然，该文给出了通过乘以某一系数的方法提高准确度的方法，是有效的。但本文推荐的方法是展开更高阶的泰勒级数，而展开的泰勒级数的阶数并不固定，而是和精度要求，以及被开方数有关。被开方数与1的接近程度可以决定泰勒级数在展开一定阶数后剩余项的大小，且可以用数学算式得出。所以，根据精度的要求，被开方数只要在0到2之间，就可以计算出应展开的阶数，然后无需乘系数，就计算出符合精度要求的平方根。

一、平方根函数的泰勒展开式

对于平方根函数 $\sqrt{x}$ ，其一阶导数为 ${d({\sqrt{x}})\over dx}={1\over 2}x^{-{1\over 2}}$ ，二阶导数为 ${d({1\over 2}x^{-{1\over 2}})\over dx}={1\over 2}\times ({-{1\over 2}})\times x^{-{3\over 2}}$ ，三阶导数为 ${d({1\over 2}\times ({-{1\over 2}})\times x^{-{3\over 2}})\over dx}={{1\over 2}\times ({-{1\over 2}})\times ({-{3\over 2}})}\times x^{-{5\over 2}}$ 。以此类推， $N$ 阶导数是 ${1\over 2}\times \prod_{c=1}^{N-1}(-{2c-1\over 2}) \times x^{-{2N-1\over 2}}$

对于泰勒展开，这里取 $x_0=1$ ，则 $x_0^{p}=1$ ，无论 $p$ 如何取值。

所以

${\sqrt{x}}=1+{x-1\over 2\times 1!}-{1\over 2}{(x-1)^2\over 2\times 2!}+{1\over 2}\times {3\over 2}{(x-1)^3\over 2\times 3!}-{1\over 2}\times{3\over 2}\times {5\over 2}{(x-1)^4\over 2\times 4!}+. . .$

由于展开式是正负交替，这里将其写成

$=1+{x-1\over 2}+\sum_{n=1}^\infty(-{\prod_{k=1}^{2n-1}(2k-1)(x-1)^{2n}\over 2^{2n}(2n)!}+ {\prod_{k=1}^{2n}(2k-1)(x-1)^{2n+1}\over 2^{2n+1}(2n+1)!})$

$={1\over 2}+{x\over 2}+\sum_{n=1}^\infty {\prod_{k=1}^{2n-1}(2k-1)(x-1)^{2n}(-2(2n+1)+(4n-1)(x-1))\over 2^{2n+1}(2n+1)!}$

$={1\over 2} + {x\over 2}+{\sum_{n=1}^\infty {\prod_{k=1}^{2n-1}(2k-1)(x-1)^{2n}(-8n+4nx-x-1)\over 2^{2n+1}(2n+1)!}}$ (1)

二、泰勒展开式收敛性的证明

式（1）中，令 $g(n) = {\prod_{k=1}^{2n-1}(2k-1)(x-1)^{2n}(-8n+4nx-x-1)\over 2^{2n+1}(2n+1)!}$

这里采用[4]中提到的比值审敛法判断式（1）的收敛性。

$\lim_{n\rightarrow \infty}{g(n+1)\over g(n)}$

$=\lim_{n\rightarrow \infty}{{\prod_{k=1}^{2(n+1)-1}(2k-1)(x-1)^{2(n+1)}(-8(n+1)+4(n+1)x-x+1)\over 2^{2(n+1)+1}(2(n+1)+1)!}\over {\prod _{k=1}^{2n-1}(2k-1)(x-1)^{2n}(-8n+4nx-x+1)\over 2^{2n+1}(2n+1)!}}$

$=\lim_{n\rightarrow \infty}{2^{2n+1}(2n+1)!\over 2^{2n+3}(2n+3)!}\times{\prod_{k=1}^{2n+1}(2k-1)(x-1)^{2n+2}(-8n+4nx+3x-7)\over \prod_{k=1}^{2n-1}(2k-1)(x-1)^{2n}(-8n+4nx-x+1)}$

$=\lim_{n\rightarrow \infty}{(2(2n)-1)(2(2n+1)-1)(x-1)^2(-8n+4nx+3x-7)\over 2^2(2n+2)(2n+3)(-8n+4nx-x+1)}$

$=\lim_{n\rightarrow \infty}{(4n-1)(4n+1)(x-1)^2(-8n+4nx+3x-7)\over 4(2n+2)(2n+3)(-8n+4nx-x+1)}$

$=\lim_{n\rightarrow \infty}{(4-{1\over n})(4+{1\over n})(x-1)^2(-8+4x+{3x\over n}-{7\over n})\over 4(2+{2\over n})(2+{3\over n})(-8+4x-{x\over n}+{1\over n})}$

$={4\times 4 \times (x-1)^2 \times (4x-8)\over 4\times 2\times 2\times (4x-8)}=(x-1)^2$

注：因为这里前提是 $x\in(0,2)$ ，所以 $4x-8$ 可以约掉。

显然，当 $x\in(0,2)$ 时， $(x-1)^2\in[0,1)$ ，故式（1）收敛。

三、泰勒展开式的阶数及剩余项的大小

（一）级数公比初步分析

首先，根据 $g(n+1)\over g(n)$ 的算式，用GeoGebra软件[5]，画出当 $x=0.5$ 和 $x=1.5$ 时 $g(n+1)\over g(n)$ 关于 $n$ 的函数图像。

从图中可看出 $g(n+1)\over g(n)$ 在 $n\geq 1$ 时单调递增，且始终大于0。鉴于其极限值存在，因此可以猜测，当 $n\geq 1$ 时， $0<{g(n+1)\over g(n)}<(x-1)^2$ 。若证实成立，则可用无穷等比数列的算式找到式（1）的剩余项的上限值。具体说明将在本章节第三小节详述。

（二）数学证明

这里用作差法证明 ${g(n+1)\over g(n)} < (x-1)^2$ 。

${g(n+1)\over g(n)} - (x-1)^2$

$={(4n-1)(4n+1)(x-1)^2(-8n+4nx+3x-7)\over 4(2n+2)(2n+3)(-8n+4nx-x+1)}-(x-1)^2$

$={(x-1)^2[(4n-1)(4n+1)(-8n+4nx+3x-7)-4(2n+2)(2n+3)(-8n+4nx-x+1)]\over 4(2n+2)(2n+3)(-8n+4nx-x+1)}$

$={(x-1)^2(-16\times 8n^3+64n^3x+48n^2x-16\times 7n^2+8n-4nx-3x+7+16\times 8n^3-64n^3x+16n^2x-16n^2+320n^2-160n^2x+40nx-40n+192n-96nx+24x-24)\over 4(2n+2)(2n+3)(-8n+4nx-x+1)}$

$={(x-1)^2[(-96n^2-60n+21)x+(192n^2+160n-17)]\over 4(2n+2)(2n+3)(-8n+4nx-x+1)}$ （2）

现在判断（2）式的正负

首先，已知 $x\in (0,2)$ ，所以 $(x-1)^2>0$ 。另外，由于 $n\geq 1$ ，所以 $4(2n+2)(2n+3)>0$ 。

对于 $(-96n^2-60n+21)x+(192n^2+160n-17)$ ，当 $n$ 确定时，该式就是 $x$ 的一次函数，在 $x\in(0,2)$ 上单调，所以只要 $x=0$ 和 $x=2$ 时该式值均大于0，该式就必然大于0。

当 $x=0$ 时，该式 $=192n^2+160n-17$ ，当 $n\geq 1$ 时其必然 $>0$

当 $x=2$ 时，该式 $=40n+25$ ，当 $n\geq 1$ 时其必然 $>0$

而对于 $-8n+4xn-x+1=(4n-1)x+1-8n$

当 $x=0$ 时，该式 $=1-8n<0$

当 $x=2$ 时，该式 $=-1<0$

所以式（2）必为负，说明 ${g(n+1)\over g(n)} < (x-1)^2$ 成立

于此同时，可用类似的方法证明 ${-8n+4nx+3x-7\over -8n+4nx-x+1}>0$ ，所以 ${g(n+1)\over g(n)}>0$

综上所述， $0<{g(n+1)\over g(n)}<(x-1)^2$

（三）剩余项上限值算式

在上一小节中，已经证明了 $0<{g(n+1)\over g(n)}<(x-1)^2$ ，所以 $\sum_{n=m}^{\infty}g(n)$ 的绝对值必然小于首项为 $g(m)$ ，公比为 $(x-1)^2$ 的无穷等比数列的绝对值。根据[6]中的无穷等比数列公式， $|\sum_{n=m}^{\infty} g(n)|<{|g(m)|\over 1-(x-1)^2}$

将 $g(m)$ 代入，可知对于式（1），当 $n$ 计算到 $m-1$ 时，剩余项的绝对值 $|o(m)|=|\sum_{n=m}^{\infty}{\prod_{k=1}^{2n-1}(2k-1)(x-1)^{2n}(-8n+4nx-x-1)\over 2^{2n+1}(2n+1)!}|<{|\prod _{k=1}^{2m-1}(2k-1)(x-1)^{2m}(-8m+4mx-x-1)|\over (1-(x-1)^2)2^{2m+1}(2m+1)!}$ （3）

因此，只要 $x\in (0,2)$ ，估算平方根时，根据所需的精度 $\epsilon$ ，只需用式（3）算出令 ${|\prod_{k=1}^{2m-1}(2k-1)(x-1)^{2m}(-8m+4mx-x-1)|\over (1-(x-1)^2)2^{2m+1}(2m+1)!}\leq \epsilon$ 的最小 $m$ 即可。

四、实验

在这一章节，用Matlab计算当 $x=0.1,0.2,. . . ,1.9$ 时通过式（1），令 $n=1,. . . ,c$ 估算出来的 $\sqrt{x}$ ，以及通过式（3）计算的 $|o(m)|$ 的上限值，因此验证通过泰勒级数算出的 $\sqrt{x}$ 估计值和实际值的误差小于（3）中计算的上限值。另外，也观察对于各 $x$ 值，通过泰勒级数展开 $c$ 阶后，余式的上限值，从而找出对于不同的 $x$ 值，需要展开多少阶后才能达到所需的精度。

（一）实验代码

首先是通过泰勒级数估算 $\sqrt{x}$ 值的程序

function res = estiSquareRoot(x, n)
    if x <= 0 || x >= 2
        res = nan;
        disp("x must be between 0 and 2, open range")
        return
    end
    if n < 0 || rem(n, 1) ~= 0 
        res = nan;
        disp("n must be a non-negative integer!")
        return
    end
    res = 1/2 + x/2;
    for m = 1:n
        ks = 1:(2 * m - 1);
        res = res + prod(2*ks-1) * (x - 1) ^ (2 * m) * (-8 * m + 4 * m * x - x - 1) / (2 ^ (2 * m + 1) * factorial(2 * m + 1));
        
    end
end

函数中的x代表 $x$ ，函数中的n代表 $c$ 。

另外，计算 $o(m)$ 的上限值

function remainSup = getRemainTaylor(x, m)
    if x <= 0 || x >= 2
        remainSup = nan;
        disp("x must be between 0 and 2, open range")
        return
    end
    if m < 0 || rem(m, 1) ~= 0 
        remainSup = nan;
        disp("n must be a non-negative integer!")
        return
    end
    ks = 1:(2*m-1);
    remainSup = prod(2 * ks - 1) * (x - 1)^(2 * m) * (-8 * m + 4 * m * x - x - 1)/((1 - (x - 1) ^ 2) * 2 ^ (2 * m + 1)*factorial(2 * m + 1));
end

函数中x代表 $x$ ，m代表 $m$ 。

然后计算对于不同的 $x$ ，不同的 $c$ ，估算的平方根值和剩余值上限。注意：对于用 $c$ 阶展开的估值，其余式是 $o(c+1)$ 。另外，也用Matlab内置的函数计算平方根，作为参考。

x_to_test = 0.1:0.1:1.9;
c_to_test = 1:10;
square_root_estimation = zeros(length(x_to_test), length(c_to_test));
square_root_remaining = zeros(length(x_to_test), length(c_to_test));
for ii_for_x = 1:length(x_to_test)
    for jj_for_c = 1:length(c_to_test)
        square_root_estimation(ii_for_x, jj_for_c) = estiSquareRoot(x_to_test(ii_for_x), c_to_test(jj_for_c));
        square_root_remaining(ii_for_x, jj_for_c) = getRemainTaylor(x_to_test(ii_for_x), c_to_test(jj_for_c) + 1);
    end
end
square_root_real = sqrt(x_to_test)';

（二）实验结果

然后观察两个结果。第一个，是估算的平方根和实际平方根的差值。

x_labels = "x=" + x_to_test;
c_labels = "c=" + c_to_test;
T = array2table(abs(square_root_real - square_root_estimation), 'VariableNames', c_labels, 'RowNames', x_labels)

结果如下：

从表格中可看出， $x$ 越远离1，精确估计平方根需要展开的泰勒阶数 $c$ 越大。

第二个，是余式的上限

T = array2table(abs(square_root_remaining), 'VariableNames',c_labels,'RowNames',x_labels)

结果如下：

比较两表，可以确认平方根的计算差距小于通过（3）计算的上限。另外，对于 $x\in[0.7,1.4]$ ， $c=1$ 时，其差距就已经小于0.001。此时其为三阶展开。对于[3]中使用的方法， $\sqrt{6}$ 如写成 $2\sqrt{1+{1\over 2}}$ ，则令 $c=1$ ，此时其根式部分的计算差不超过0.0021，整体差距不超过0.0042。同样，若写成 $3\sqrt{1-{1\over 3}}$ ，则根式部分计算差不超过0.0015，整体差距不超过0.003。当然，若 $x$ 离1更远，需展开的阶数会更大，但只要 $x\in(0,2)$ ，都可通过式（3）算出需展开的阶数。