数值分析实验 - Hamming 级数求和

数值分析实验 1 - Hamming 级数求和

实验要求

\[ \psi(x)=\sum_{k=1}^{\infty}\frac{1}{k(k+x)}\tag{1} \]

\(x=0.0,~0.1,~0.2~,...,300.0\), 计算在这 \(3001\) 个值下所有的 \(\psi(x)\), 并将误差控制在 \(1.0\rm{e-}10\) 范围

分析

直接计算的耗时往往过长 (见 直接计算的误差及复杂度分析), 故需要进行优化

本文给出了一种将绝大多数数据通过线性递推计算,从而以极快的运算速度完成运算的算法

计算公式

注意到

  • \[ \psi(0)=\sum_{k=1}^{\infty}\frac{1}{k^2}=\frac{\pi^2}{6}\tag{2} \]

  • \(x\ne0\)

    \[ \psi(x)=\sum_{k=1}^{\infty}\frac{1}{k(k+x)}=\frac{1}{x}\sum_{k=1}^{\infty}\left(\frac{1}{k}-\frac{1}{k+x}\right)\tag{3} \]

    注意到,当 \(x\geqslant1\)

    \[ \begin{aligned} x~\psi(x)&=\sum_{k=1}^{\infty}\frac{1}{k}-\sum_{k=1}^{\infty}\frac{1}{k+x}\\ &=\sum_{k=1}^{\infty}\frac{1}{k}-\sum_{k=1}^{\infty}\frac{1}{k+x-1}+\frac{1}{x}\\ &=(x-1)~\psi(x-1)+\frac{1}{x} \end{aligned} \]

    则有

    \[ \psi(x)=\frac{x-1}{x}\psi(x-1)+\frac{1}{x^2}\tag{4} \]

对于 \(0<x<1\) 的点,我们这样考虑

注意到

\[ \sum_{k=a+1}^{\infty}\frac{1}{k^m}<\int_a^{\infty}{\mathrm{d}k\over k^m}\tag{5} \]

我们需要取满足

\[ \int_a^{\infty}{\mathrm{d}k\over k^m}\leqslant\epsilon=10^{-10}\tag{6} \]

的最小 \(a\) 值,且让其尽可能小 (至少小于 \(10^7\))

容易解得

\[ a\geqslant(\epsilon(m-1))^\frac{1}{1-m}=\sqrt[m-1]{10^{10}\over m-1},\qquad m>1\tag{7} \]

\(m=3\) 时,\(a\geqslant 7.0711\times10^4\)

这启发我们应尝试将级数 \(\psi(x)\) 分母项的次数升高到 \(3\)

注意到

  • \[ \psi(1)=\sum_{k=1}^{\infty}\frac{1}{k(k+1)}=1\tag{8} \]

我们设级数

\[ \psi_1(x)=\frac{\psi(x)-\psi(1)}{1-x}=\sum_{k=1}^{\infty}\frac{1}{k(k+1)(k+x)}\tag{9} \]

则可将

\[ (1-x)\sum_{k=1}^{a_3}\frac{1}{k(k+1)(k+x)}+1\tag{10} \]

视为 \(\psi(x)\) 的近似值,其中 \(a_3\geqslant 70711\), 下面代码中取 \(71000\)

所以最终的计算公式如下

\[ \psi(x)=\begin{cases} \frac{\pi^2}{6},&x=0\\ (1-x)\sum_{k=1}^{a_3}\frac{1}{k(k+1)(k+x)}+1,&0<x<1\\ \frac{x-1}{x}\psi(x-1)+\frac{1}{x^2},&x\geqslant 1 \end{cases}\tag{11} \]

时间复杂度为 \(\Theta(a_3+n)\)

程序设计

C++ 程序

Show code

main.cppview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
/*
* @Author: Tifa
* @LastEditTime: 2021-04-16 20:30:47
* @Description: Exp.1
*/
// -std=c++14
#include <cmath>
#include <cstdio>

using num_t = double;

// data
const int N = 3001;
const num_t pi = acos(-1.0);
// iteration
const int A_3 = 71000;
// step
const num_t STEP = 1e-1;
// inverse of step
const num_t INV_STEP = 1.0 / STEP;

// result
num_t psi[N] = {pi * pi / 6};

int main() {
// 0 < x < 1
for (int i = 1; i < INV_STEP; ++i) {
num_t psi_now = 0;
for (int j = 1; j <= A_3; ++j)
psi_now += 1.0 / (j * (j + 1) * (j + STEP * i));
psi[i] = psi_now * (1 - STEP * i) + 1;
}
// x >= 1
for (int i = INV_STEP; i < N; ++i)
psi[i] = (1 - INV_STEP / i) * psi[(size_t)(i - INV_STEP)] +
INV_STEP * INV_STEP / (i * i);
for (int i = 0; i < N; ++i) printf("%6.2f %16.12f\n", STEP * i, psi[i]);
return 0;
}

MATLAB 程序

主程序

Show code

main.mview raw
1
2
3
4
5
6
7
% Exp.1

% @Author: Tifa
% @LastEditTime: 2021-04-16 20:30:47

A = [0:0.1:300; calc_fast()'];
sprintf('%6.2f %16.12f\n', A)

计算函数

Show code

calc_fast.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
function ret = calc_fast()
% Syntax: ret = calc_fast()
%
% The function about calculating Hamming's series

% @Author: Tifa
% @LastEditTime: 2021-04-16 20:30:47

% check results when this variable is true
check = true;

ret = zeros(3001, 1);
ret(1) = pi^2/6;

syms k s

for x = 0.1:0.1:0.9
s = 1 / (k * (k + 1) * (k + x));
ret(floor(x * 10) + 1) = (1 - x) * vpa(symsum(s, k, 1, 71000)) + 1;
end

for x = 1:0.1:300
ret(floor(x * 10) + 1) = (x - 1) * ret(floor(x * 10) - 9) / x + 1 / x / x;
end

if check

for x = 0:0.1:300

if error_judge(x, ret(floor(x * 10) + 1))
error('Error exceed!')
end

end

end

end

误差检验函数

Show code

error_judge.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
function output = error_judge(x, input_sum)
% Syntax: output = error_judge(x, input_sum)
%
% return true if the input value's error exceeded, otherwise return false

% @Author: Tifa
% @LastEditTime: 2021-04-16 20:30:47

syms k s
s = 1 / (k * (k + x));
accurate_sum = symsum(s, k, 1, inf);

output = abs(input_sum - vpa(accurate_sum)) > 1e-10;
end

结果讨论和分析

直接计算的误差及复杂度分析

\[ \psi^*(x)=\sum_{k=1}^a\frac{1}{k(k+x)}\tag{12} \]

注意到

  • \[ \sum_{k=a+1}^{\infty}\frac{1}{k(k+x)}<\int_a^{\infty}{\mathrm{d}k\over k(k+x)}=\frac{\ln(a+x)-\ln a}{x}\tag{13} \]
  • \[ \frac{\ln(a+x)-\ln a}{x}\leqslant\epsilon\implies a\geqslant{x\over e^{x\epsilon}-1}\tag{14} \]
  • \(x\ll\frac{1}{\epsilon}\) 时,\[ {x\over e^{x\epsilon}-1}\approx\frac{1}{\epsilon}\tag{15} \]

本实验中,\(x_{max}=300.0,~\epsilon=1.0\rm{e-}10\)

\(a\) 大致取 \(\frac{1}{\epsilon}\) 即可保证 \(\psi(x)\) 的绝对误差 \(|E(\psi(x))|\leqslant\epsilon\), 实验结果也符合这一估计

则直接计算的时间复杂度为

\[ \Theta\left(\sum_{x=1}^n\left(0.1x+\frac{1}{\epsilon}\right)\right)=\Theta\left(\frac{n}{\epsilon}+n^2\right) \]

其中

  • \(n\)\(x\) 样本点的个数,本实验中即为 \(3001\)
  • \(\epsilon\) 为误差范围,本实验中即为 \(1.0\rm{e-}10\)

可以发现,由于 \(\epsilon\) 很小,所以直接计算的耗时会很长。在每秒执行 \(1\rm{e}8\) 次运算的计算机上预计花费一天多的时间才能完成计算

进一步的优化方案

  • 通过提高 \(m\) 来减小 \(a\) 我们知道,在满足误差要求的条件下,\(a\)\(m\) 满足如下关系 (即式 \((7)\)):

    \[ a\geqslant\sqrt[1-m]{\epsilon(m-1)}=\sqrt[m-1]{10^{10}\over m-1},\qquad m>1 \]

    容易得知,\(f(x)=\sqrt[1-m]{\epsilon(m-1)},~(m>1)\) 是严格递减的凸函数,且有下表

    \(m\) \(\lceil f(m)\rceil\) \(m\lceil f(m)\rceil\)
    \(2\) \(10^{10}\) \(2\times 10^{10}\)
    \(3\) \(70711\) \(212133\)
    \(4\) \(1494\) \(5976\)
    \(5\) \(224\) \(1120\)
    \(6\) \(73\) \(438\)
    \(7\) \(35\) \(245\)
    \(8\) \(21\) \(168\)

    由于计算每个数的时候都会做至少 \(m\lceil f(m)\rceil\) 次乘法,且考虑累计误差的影响,我们只考虑 \(m=4\) 的情况

    \[ \psi_2(x)=\frac{\psi_1(x)-\psi_1(2)}{2-x}\tag{16} \]

    由式 \((9)\), 有

    \[ \psi_2(x)=\sum_{k=1}^{\infty}\frac{1}{k(k+1)(k+2)(k+x)}\tag{17} \]

    则可将

    \[ (1-x)\left((2-x)\sum_{k=1}^{a_4}\frac{1}{k(k+1)(k+2)(k+x)}+\frac{1}{4}\right)+1\tag{18} \]

    视为 \(\psi(x)\) 的近似值,其中 \(a_4\geqslant 1494\), 推荐取 \(1500\)