*[[Programable Sourceの使い方]] [#a59dbce3]

可視化ソフトウェア[[ParaView:http://www.paraview.org]]では、"Programmable Source"を使用して、データを作成することができます。

ここでは、"Programmable Source"でImageデータ(Uniform Rectilinear Grid)を作成してみます。データの各点には、Riemann Zeta関数の絶対値を与えます。

**Programmable Sourceの場所 [#f0a44582]

Programmable Sourceはメニューバーの
  Sources >> Programmable Source
にあります。

#ref(menue_programmable_source.png,center,nowrap,40%,programmable source);

**Script (RequestInformation)の記述 [#gbe5d860]

まず、300x500x1ピクセルの2次元画像を作成します。

Propertiesタブにある"Output Data Set Type"に"vtkImageData"を設定し、"Script (RequestInformation)"に以下を記述します。

#code(python){{
from paraview import util
util.SetOutputWholeExtent(self, [0, 299, 0, 499, 0, 0])
}}

 
#ref(RequestInformation.png,center,nowrap,40%,RequestInformation);

**Scriptの記述 [#f0bef0f2]

Propertiesタブにある"Script"には以下を記述します。

これは、画像データ上の点でリーマン・ゼータ関数を求めて、その絶対値をその点に設定します。Riemann Zeta関数は、

\[
\zeta(s) = 1 + \frac{1}{2^s} + \frac{1}{3^s} + \frac{1}{4^s} + \cdots
\]

で定義される無限級数の複素関数です。

詳しくは、[[リーマンのゼータ関数で遊び倒そう (Ruby編):http://tsujimotter.hatenablog.com/entry/riemann-zeta-function]]を参照のこと。

なお、メッシュの次元(300, 500, 1)、原点(0.0001, 0.0001, 0)、間隔(0.1, 0.1, 1)はコード内の以下で設定しています。

  out.SetDimensions(300, 500, 1)
  out.SetOrigin(0.0001, 0.0001, 0)
  out.SetSpacing(0.1, 0.1, 1)

#code(python){{
import math

lowerThreshold = 1.0e-6
upperBound = 1.0e+4
infinity = 300

def binom(n,k):
  m = 1
  if n < 2 * k:
    k = n - k
  for i in range(1, k + 1):
    m = m * (n - i + 1) / i
  return m

def zeta(s):

  outerSum= 0.0
  prev = 1000000000

  for m in range(1, infinity+1):

    innerSum = 0.0
    for j in range(1, m+1):
      c1 = 1 if (j-1)%2 == 0 else -1
      c2 = binom(m-1, j-1)
      c3 = j**(-s)
      innerSum += c1*c2*c3

    try:
      innerSum = innerSum * (2**(-m)) / (1-2**(1-s))
    except:
      return 0+0j

    outerSum += innerSum

    if abs(prev - innerSum) < lowerThreshold:
      break

    if abs(outerSum) > upperBound:
      break

    prev = innerSum

  return outerSum


out = self.GetOutput()
out.SetDimensions(300, 500, 1)
out.SetOrigin(0.0001, 0.0001, 0)
out.SetSpacing(0.1, 0.1, 1)

zetaArray = vtk.vtkDoubleArray()
zetaArray.SetName('zeta')
zetaArray.SetNumberOfComponents(1)
zetaArray.SetNumberOfTuples(out.GetNumberOfPoints())

dims = out.GetDimensions()

for z in range(0, dims[2]):
  for y in range(0, dims[1]):
    for x in range(0, dims[0]):
      pntId = out.ComputePointId([x, y, z])
      point = out.GetPoint(pntId)
      s = point[0] + point[1]*1j
      val = abs(zeta(s))
      zetaArray.SetTuple1(pntId, val)

out.GetPointData().AddArray(zetaArray)
}}

これは、画像データ上の点でリーマン・ゼータ関数を求めて、その絶対値をその点に設定します。Riemann Zeta関数は、

\[
\zeta(s) = 1 + \frac{1}{2^s} + \frac{1}{3^s} + \frac{1}{4^s} + \cdots
\]

で定義される無限級数の複素関数です。

詳しくは、[[リーマンのゼータ関数で遊び倒そう (Ruby編):http://tsujimotter.hatenablog.com/entry/riemann-zeta-function]]を参照のこと。

なお、メッシュの次元(300, 500, 1)、原点(0.0001, 0.0001, 0)、間隔(0.1, 0.1, 1)はコード内の以下で設定しています。

  out.SetDimensions(300, 500, 1)
  out.SetOrigin(0.0001, 0.0001, 0)
  out.SetSpacing(0.1, 0.1, 1)

 
#ref(script.png,center,nowrap,40%,script);

上記の設定が完了しましたら、"Apply"します。(Riemann Zeta関数の計算は時間がかかりますので気長に待ちましょう)。

**描画の設定 [#r5185f5d]

Riemann Zeta関数の計算が終了し、画像データが作成されましたら、表示を変数"zeta"にしコンター図を表示しましょう。

ここで、Riemann Zeta関数は、s=1が極となり、その周辺で絶対値が非常に大きいため、コンター表示する場合は、値の範囲を"Rescale"する必要があります。

また、形状に座標軸(実部(real)と虚部(imaginary))を表示するとよいでしょう。

以下の図のように、値の範囲を[0, 4]にリスケールし、続けて、Axes Gridsを有効にしてから、"Edit"ボタンから軸の名前を(x軸のタイトル->real, y軸のタイトル->imaginary)に変更します。Axes GridsのEditでは、フォントの種類や大きさも変更できますので適宜調節すると良いでしょう。

#ref(rescale_and_axis_grid.png,center,nowrap,40%,rescale and enable axis_grid);

以上で、Programable Sourceを使用して、画像データを作成し、その各点にRiemann Zeta関数の絶対値を与え、最期に、そのコンターを見ることができました。

#ref(riemann_zeta_contours.png,center,nowrap,40%,riemann-zeta contours);

**Riemann予想 [#fb775def]

Riemann予想によるとRiemann Zeta関数の零点

\[
\zeta(s) = 1 + \frac{1}{2^s} + \frac{1}{3^s} + \frac{1}{4^s} + \cdots = 0を満たすsの値
\]

の実部は全て1/2となります。

Slice機能で、実部が1/2の位置で切断して、表示してみると確かにその線上に、零点(青色)が並んでいます。

#ref(Riemann_zero_point.png,center,nowrap,40%,real part=0.5);

また、Warp by Scalarフィルターなどを使うと、零点の位置が見易いかもしれません。

#ref(warp_by_scalar.png,center,nowrap,30%,warp by scalar);

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS