找回密码
 欢迎注册
查看: 17520|回复: 3

[提问] Project Euler Problem 368 modified harmonic series

[复制链接]
发表于 2012-1-23 22:08:01 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?欢迎注册

×
Problem 368
22 January 2012



The harmonic series 1 + 1 /2
+ 1 /3
+ 1 /4
+ ... is well known to be divergent.  


If we however omit from this series every term where the denominator has a 9 in it, the series remarkably enough converges to approximately 22.9206766193.
This modified harmonic series is called the Kempner series.

Let us now consider another modified harmonic series by omitting from the harmonic series every term where the denominator has 3 or more equal consecutive digits. One can verify that out of the first 1200 terms of the harmonic series, only 20 terms will be omitted.
These 20 omitted terms are:

1/111
, 1/222
, 1 /333
, 1 /444
, 1 /555
, 1 /666
, 1 /777
, 1 /888
, 1 /999
, 1 /1000
, 1 /1110
, 1/1111
, 1 /1112
, 1 /1113
, 1 /1114
, 1 /1115
, 1 /1116
, 1 /1117
, 1 /1118
and  1/1119
.

This series converges as well.

Find the value the series converges to.
Give your answer rounded to 10 digits behind the decimal point.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2012-1-27 11:49:09 | 显示全部楼层
这个应该可以用PolyGamma函数的展开式帮助计算:
http://bbs.emath.ac.cn/viewthrea ... fromuid=20#pid28839
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2012-1-28 21:27:00 | 显示全部楼层
网上找到一个代码,谁能解释下

  1. function [ s,error ] = Problem42
  2. % Thomas Schmelzer, November 2008
  3. % computes the limit of the Kempner series \sum_{n=1}^\infty 1/n
  4. % where n does not contain the substring "42" --- a problem that has been
  5. % proposed by Folkmar Bornemann. For details see
  6. % Schmelzer, Thomas; Robert Baillie (June朖uly 2008). "Summing a Curious,
  7. % Slowly Convergent Series". American Mathematical Monthly 115 (6): 525?40
  8. % exact solution: 228.44630 41592 30813 25414 80861 26250 58957 81629 ...
  9. %% given the recurrence matrix T, compute f,A and B.
  10. T = [1 1 1 1 2 1 1 1 1 1; 1 1 0 1 2 1 1 1 1 1];
  11. [I,J] = find(T); f = zeros(2,2,10); A = zeros(2,2);
  12. for s = 1:length(I)
  13.     f(T(I(s),J(s)),I(s),J(s))=1;
  14. end
  15. for s = 0:9
  16.     A = A + f(:,:,s+1)/10;
  17. end
  18. I = eye(size(A)); B = inv(I-A)-I;

  19. %% define the S_1 and S_2 explicitly for integers with up to 5 digits
  20. S = cell(5,2); S{1,1} = [1,2,3,5,6,7,8,9]; S{1,2} = 4;
  21. for i = 2:size(S,1)
  22.     for m = 0:1:9    % possible digits to attend
  23.         [J,L] = find(f(:,:,m+1));
  24.         for s = 1:length(J)
  25.             S{i,J(s)} = [S{i,J(s)},10*S{i-1,L(s)}+m];
  26.         end
  27.     end
  28. end

  29. %% define the numerical parameters
  30. K = 20;  % extrapolation   
  31. P = 30;  % power cutoff
  32. Ignore = -22;

  33. Psi = zeros(K,2,P);
  34. % compute the sums Psi_{i} explicitly for up to 5 digits
  35. for i = 1:size(S,1)
  36.     for k = 1:P
  37.         Psi(i,1,k) = sum(S{i,1}.^(-k));
  38.         Psi(i,2,k) = sum(S{i,2}.^(-k));
  39.     end
  40. end

  41. warning('off','MATLAB:log:logOfZero');
  42. for i = size(S,1)+1:1:K
  43.     for k = 1:P
  44.         if (Psi(i-1,1,k)>0)
  45.             for m = 0:9
  46.                 [J,L] = find(f(:,:,m+1));
  47.                 for s = 1:length(J)
  48.                     for w = 0:P-k
  49.                        if log10(abs(aCoeff(k,w,m))) + log10(abs(Psi(i-1,L(s),k+w))) > Ignore;
  50.                            Psi(i,J(s),k) = Psi(i,J(s),k)+aCoeff(k,w,m)*Psi(i-1,L(s),k+w);
  51.                        end        
  52.                     end
  53.                 end
  54.             end   
  55.         end
  56.     end
  57. end
  58. warning('on','MATLAB:log:logOfZero');
  59. % Extrapolation
  60. s = sum(B*Psi(K,:,1)') + sum(sum(Psi(:,:,1)));
  61. error = 228.44630415923081325414 - s;

  62. function [ a ] = aCoeff( k,w,m )
  63. if (m+w==0)
  64.     mw = 1;
  65. else
  66.     mw = m^w;
  67. end
  68. a = 10^(-k-w)*mw*(-1)^w*nchoosek(k+w-1,w);
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2012-1-28 21:30:33 | 显示全部楼层
还有这个,但不知怎么改造成368的答案
kempnerSums.m (56.28 KB, 下载次数: 3)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-5-4 22:53 , Processed in 0.050344 second(s), 19 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表