评估总和

计算科学 表现 正则 C 特殊功能
2021-12-23 18:24:52

我想评估总和

k=1(i+12)kkα
在哪里i=1α[34,1]具有 8 位精度。

如果我愿意每次计算最多花费一秒钟的 CPU 时间,我应该如何进行?

假设我想要一个非常快速的评估方案并且我愿意花费大量的工作。我应该如何进行?

任何线索、线索或建议将不胜感激。

4个回答

你试过什么了?这个完全幼稚的实现设法在我的笔记本电脑上在 2.5 秒内计算出 7 个(可能是 7.5 个)数字:

#include <iostream>
#include <complex>
#include <cmath>
#include <iomanip>

int main ()
{
  const double alpha = 1;
  std::cout.precision(16);

  std::complex<double> sum = 0;
  for (unsigned int k=1; k<10000000; ++k)
    {
      sum += std::pow(std::complex<double>(1,1)/std::sqrt(2.), k)
         *
         std::pow(k, -alpha);

      if (k % 1000000 == 0)
    std::cout << k << ' ' << sum << std::endl;
    }
}

结果是:

g/x> c++ -O3 x.cc
g/x> time ./a.out 
1000000 (0.2674004983680959,1.178096037989336)
2000000 (0.2674002483693692,1.17809664154278)
3000000 (0.2674001650362904,1.178096842727255)
4000000 (0.267400123369689,1.178096943319499)
5000000 (0.2674000983697316,1.178097003674927)
6000000 (0.2674000817030793,1.178097043911768)
7000000 (0.2674000697983399,1.178097072652377)
8000000 (0.2674000608697807,1.17809709420789)
9000000 (0.2674000539253223,1.178097110973249)

real    0m2.553s
user    0m2.548s
sys     0m0.001s

我确信这可以大大改善,例如通过使用以下事实(1+i)k只能获得 8 个不同的值,可以制成表格,即2k可以很容易地计算出来,等等。我毫不怀疑,通过 20 分钟的工作,我可以在一秒钟内完成 8 位数的计算。对于较小的值可能更难做α,然而,收敛速度较慢。

由于该问题有一个 Fortran 标签,因此这是 Wolfgang 在 Fortran 中的解决方案:

implicit none
integer, parameter :: dp = kind(0.d0)
complex(dp), parameter :: i_ = (0, 1)

real(dp) :: alpha = 1
complex(dp) :: s = 0
integer :: k
do k = 1, 10000000
    s = s + ((i_+1)/sqrt(2._dp))**k * k**(-alpha)
    if (modulo(k, 1000000) == 0) print *, k, s
end do
end

结果是:

$ gfortran -O3 x.f90 
$ time ./a.out 
     1000000 ( 0.26740049836809593     ,  1.1780960379893362     )
     2000000 ( 0.26740024836936921     ,  1.1780966415427796     )
     3000000 ( 0.26740016503629038     ,  1.1780968427272547     )
     4000000 ( 0.26740012336968905     ,  1.1780969433194985     )
     5000000 ( 0.26740009836973161     ,  1.1780970036749265     )
     6000000 ( 0.26740008170307927     ,  1.1780970439117682     )
     7000000 ( 0.26740006979833991     ,  1.1780970726523770     )
     8000000 ( 0.26740006086978074     ,  1.1780970942078899     )
     9000000 ( 0.26740005392532235     ,  1.1780971109732488     )
    10000000 ( 0.26740004836978204     ,  1.1780971243856078     )

real    0m1.477s
user    0m1.472s
sys 0m0.000s

为了比较,在我的计算机上,C++ 代码采用:

$ g++ -O3 x.cc
$ time ./a.out 
1000000 (0.2674004983680959,1.178096037989336)
2000000 (0.2674002483693692,1.17809664154278)
3000000 (0.2674001650362904,1.178096842727255)
4000000 (0.267400123369689,1.178096943319499)
5000000 (0.2674000983697316,1.178097003674927)
6000000 (0.2674000817030793,1.178097043911768)
7000000 (0.2674000697983399,1.178097072652377)
8000000 (0.2674000608697807,1.17809709420789)
9000000 (0.2674000539253223,1.178097110973249)

real    0m2.123s
user    0m2.116s
sys 0m0.004s

我使用 gcc 4.6.3。

mpmath 有这个

>>> 导入 mpmath
>>> mpmath.polylog(1, (1+1j)/mpmath.sqrt(2))
mpc(真实='0.26739999836978523',图像='1.1780972450961724')
>>> 导入 numpy
>>> for x in numpy.linspace(0.75, 1.0, 11): print mpmath.polylog(x, mpmath.sqrt(1j))
...
(0.13205298513153 + 1.22168160828646j)
(0.147070653672424 + 1.21783198879211j)
(0.161744021745461 + 1.21384547018142j)
(0.17607856182228 + 1.20973163332567j)
(0.190079791459314 + 1.20549967730407j)
(0.2037532636373 + 1.20115843051452j)
(0.217104557681381 + 1.19671636161604j)
(0.230139270736483 + 1.19218159029609j)
(0.242863009773464 + 1.18756189785626j)
(0.25528138410235 + 1.18286473761121j)
(0.267399998369785 + 1.17809724509617j)
>>>

如果您绘制该系列的部分总和,您将得到以下结果:

在此处输入图像描述

这是真正的部分α=0.91对于 10 到 100 项的部分总和。显然,这里有一个模式,虚部也遵循类似的模式。似乎您可以取 8 个连续的部分和并对它们进行平均以获得更好的近似值。或者至少您可以从这些部分和中获得渐近数值行为。

更具分析性,您可以将总和分解为 8 个不同单位根上的单独总和:

j=18(i+12)jk=0(8k+j)α

每个内部和都是一个 Hurwitz zeta 函数,其中有特殊的评估例程(如GSL中)。以这种方式绘制总和作为α

在此处输入图像描述

蓝线是实部,红线是虚部。这似乎表现得很好。