Visualizing the primes with the Riemann Zeta function

Question

I am trying to plot the identity seen here, namely that if we define: $$\psi _{0}(x)={\frac 12}\left(\sum _{{n\leq x}}\Lambda (n)+\sum _{{n<x}}\Lambda (n)\right)$$ Then, it equal to the following, where the sum is taken over the nontrivial zeta zeroes. $$\psi _{0}(x)=x-\sum _{{\rho }}{\frac {x^{{\rho }}}{\rho }}-{\frac {\zeta '(0)}{\zeta (0)}}-{\frac {1}{2}}\log(1-x^{{-2}})$$ The linked website show a very nice graphic where $n$ is the number of zeta zeroes used to approximate the function:

However, I am having trouble reproducing this. My goal is just to plot this for one value of $n$. So, my exact definition of the function in question seems to be working fine:

    Ches0[x_] := 
 1/2 (Sum[MangoldtLambda[n], {n, 1, x}] + 
    Sum[MangoldtLambda[n], {n, 1, x - 1}])

However, my implementation of the approximate identity, while it gets the values at the jumps right, seems to be very off in most other respects.

Ches0Appr[x_] := 
 x - Log[2*Pi] - Sum[x^ZetaZero[i]/ZetaZero[i], {i, 1, 100}] - 
  1/2*Log[1 - 1/x^2]

I'm not exactly sure what's going wrong here, so my question is, how can I plot the second formula in my post?

Answer 1

As others have implicitly mentioned, your implementation of the second formula omits the non-trivial zeros with negative imaginary part. When the sum is taken over all zeros, then the imaginary part of the sum cancels, leaving just a real part. For example,

Simplify[ComplexExpand[Re[x^(1/2 + I*t)/(1/2 + I*t) +
                          x^(1/2 - I*t)/(1/2 - I*t)]],
         Assumptions -> x >= 1]

The result involves only the real values t and x.

(4 Sqrt[x] (Cos[t Log[x]] + 2 t Sin[t Log[x]]))/(1 + 4 t^2)

If the Psi[x,n] function is defined as follows, the calculations using the code from @ubpdqn are much faster.

ChebyshevPsiZeta[x_Integer, n_Integer] :=
   Block[{t = N[Im[ZetaZero[Range[n]]]], tlogx},
      tlogx = t*Log[x];
      x - Log[2 Pi] - Log[1 - 1/x^2]/2 - 
      4 Sqrt[x] Total[ (Cos[tlogx] + 2 t Sin[tlogx])/(1 + 4 t^2)]]

Define vis2 which includes ChebyshevPsiZeta[x,n] instead of f[x,n] in vis from @ubpdqn.

vis2[n_] := 
   Plot[{g[x], ChebyshevPsiZeta[x, n]}, {x, 1, 100},
        Exclusions -> None,
        PlotStyle -> {Green, Red}, Background -> Black, Frame -> True, 
        PlotLabel -> Style[Row[{"n=", n}], 20, White], ImageSize -> 400, 
        PlotRange -> {0, 100}]

Using Map[vis2,{1,2,3,4,5,6,7,8,9,10,20,30,40,50}] was about 25 times faster than with vis. Using ParallelMap was about 90 times faster.

Also note that the following definition of the von Mangoldt function by Stefan is faster than the built-in MangoldtLambda.

vonMangoldtLambda[x_Integer] := If[Length[#]==1, Log[#[[1,1]]], 0]&[FactorInteger[x]]
SetAttributes[vonMangoldtLambda, Listable]

Answer 2

EDIT

Based on useful J.M.'s comments: This is not efficient. However, it achieves result (I think):

f[x_, n_] := Module[{zero = Flatten[N[Table[ZetaZero[{j, -j}], {j, n}]]]}, 
x - Log[2 Pi] - Total[x^zero/zero] - Log[1 - 1/x^2]/2]
g[x_] := Sum[MangoldtLambda[j], {j, 1, Floor[x]}]

vis[n_] := 
 Plot[{g[x], f[x, n]}, {x, 1, 100}, Exclusions -> None, 
  PlotStyle -> {Green, Red}, Background -> Black, Frame -> True, 
  PlotLabel -> Style[Row[{"n=", n}], 20, White], ImageSize -> 400, 
  PlotRange -> {0, 100}]

After a 4 km (40 minute) walk the following allowed export of gif:

animation = 
  vis[#] & /@ {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 
    80, 90, 100, 150, 200, 250};