Как построить спектр дзета-нуля Римана с преобразованием Фурье в Mathematica? - PullRequest
Новая
       16

Как построить спектр дзета-нуля Римана с преобразованием Фурье в Mathematica?

12 голосов
/ 20 января 2012

В работе «Гипотеза Римана» Дж. Брайана Конри на рис. 6 приведен график преобразования Фурье слагаемого ошибки в теореме о простых числах.См. График слева на изображении ниже:

Plots from Conrey's paper on the Riemann hypothesis

В блоге под названием Простые числа из тонкого воздуха , написанном Крисом Кингом, есть MatlabПрограмма, которая строит спектр.Смотрите сюжет справа в начале поста.Возможен перевод в Mathematica:

Mathematica: 

 scale = 10^6;
 start = 1;
 fin = 50;
 its = 490;
 xres = 600;
 y = N[Accumulate[Table[MangoldtLambda[i], {i, 1, scale}]], 10];
 x = scale;
 a = 1;
 myspan = 800;
 xres = 4000;
 xx = N[Range[a, myspan, (myspan - a)/(xres - 1)]];
 stpval = 10^4;
 F = Range[1, xres]*0;

For[t = 1, t <= xres, t++,
 For[yy=0, yy<=Log[x], yy+=1/stpval,
 F[[t]] =
 F[[t]] +
 Sin[t*myspan/xres*yy]*(y[[Floor[Exp[yy]]]] - Exp[yy])/Exp[yy/2];
 ]
 ]
 F = F/Log[x];
 ListLinePlot[F]

Однако, как я понимаю, это матричная формулировка преобразования синуса Фурье, и поэтому его вычисление очень дорого.Я НЕ рекомендую запускать его, потому что он уже один раз сломал мой компьютер.

Есть ли способ в Mathematica, использующий быстрое преобразование Фурье, построить график с пиками при значениях x, равных мнимой части дзета-нулей Римана?

Я пробовал команды FourierDST и Fourier безуспешно.Кажется, проблема в том, что переменная yy в коде включена в Sin[t*myspan/xres*yy] и (y[[Floor[Exp[yy]]]] - Exp[yy])/Exp[yy/2].

РЕДАКТИРОВАТЬ: 20.1.2012, я изменил строку:

For[yy = 0, yy <= Log[x], 1/stpval++,

в следующее:

For[yy = 0, yy/stpval <= Log[x], yy++,

РЕДАКТИРОВАТЬ: 22.1.2012, Из комментария Хайке, изменено:

For[yy = 0, yy/stpval <= Log[x], yy++,

в:

For[yy=0, yy<=Log[x], yy+=1/stpval,

1 Ответ

11 голосов
/ 23 января 2012

А как насчет этого?Я слегка переписал синусоидальное преобразование, используя идентификатор Exp[a Log[x]]==x^a

Clear[f]
scale = 1000000;
f = ConstantArray[0, scale];
f[[1]] = N@MangoldtLambda[1];
Monitor[Do[f[[i]] = N@MangoldtLambda[i] + f[[i - 1]], {i, 2, scale}], i]

xres = .002;
xlist = Exp[Range[0, Log[scale], xres]];
tmax = 60;
tres = .015;
Monitor[errList = Table[(xlist^(-1/2 + I t).(f[[Floor[xlist]]] - xlist)), 
  {t, Range[0, 60, tres]}];, t]

ListLinePlot[Im[errList]/Length[xlist], DataRange -> {0, 60}, 
  PlotRange -> {-.09, .02}, Frame -> True, Axes -> False]

, который производит

Mathematica graphics

...