「数学」カテゴリーアーカイブ

プログラミングと数値計算, 数学

二重指数関数型数値積分

1件のコメント

二重指数関数型数値積分公式(DE公式またはtanh-sinh公式)とは、変数変換によって被積分関数を別の関数に変換し積分する方法のことです。
計算アイデアは

  1. 変数変換によって台形則が良く働く被積分関数に変換する
  2. 台形則によって数値積分する

という考えに立脚しています。端点に特異点が存在する場合に最適、という特徴を持ちます。
1974年に、高橋 秀俊と森 正武によって提案されました[1]。

二重指数関数型数値積分とは、関数\(f(x)\)の\([-1,1]\)に渡る積分を

\(
\begin{align}
\int_{-1}^{1} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\tanh\left(\frac{\pi}{2}\sinh(ih)\right), \\
w_i&=\frac{h\frac{\pi}{2}\cosh(ih)}{\cosh^2\left(\frac{\pi}{2}\sinh(ih)\right)}
\end{align}
\)

として近似します。\(N_-, N_+\)は離散化誤差と打ち切り誤差が等しくなるように決められます。
\(N_-, N_+\)はプログラム上では、\(x_i\)がコンピュータの扱える桁数を超えず、\(w_i\)がアンダーフローを起こさない範囲で決められます。


[adsense1]

計算アイデア

DE公式は変数変換型の数値積分公式と呼ばれます。
DE公式の計算は、

  1. 変数変換によって台形則が良く働く被積分関数に変換する
  2. 台形則によって数値積分する

という考えが元になっています。

大きな特徴として、

・端点特異性に強い
・応用範囲が広い

ことが挙げられます。
例えば、端点で発散してしまう積分
\(
\displaystyle \int_0^1 \frac{1}{\sqrt{x}} dx
\)
や、厄介な積分
\(
\displaystyle \int_0^1 \sin\left(\frac{1}{\sqrt{x}}\right)/\sqrt{x} dx
\)
の計算も少ない分点数で実行できます(後者の厄介な積分はDE公式を用いて4桁程度一致します。しかし、ほかの補間型の積分公式では1桁合うかどうかでしょう)。

なぜ変数変換でうまく計算できるのか、これを理解するにはまず、台形則が良く働く場合とは何かを知らなければなりません。

台形則が良く働く関数とは、関数の端点での微分の値が近い事です。
台形則の刻み幅を\(h\), 区間\([a,b]\)を台形則によって関数\(f(x)\)を数値積分する時、本来の積分値と台形則で求めた値の誤差を表す第1項は、
\(
\displaystyle \frac{h^2}{12}\{f'(b)-f'(a)\}
\)
と表されます。刻み幅を小さくすれば誤差が小さくなるのは直観的に分かりますが、誤差を小さくできる要因は\(\{f'(b)-f'(a)\}\)にもあります。これは、
・両端で微分の値が等しい(周期関数の1周期に渡る積分、\(a,b\)に近づくつれて関数が定数に近づいていく積分)
を考えるとき、高精度の結果を与えると言っています。

上記条件が満たされるとき、他の高次の誤差項\(f^{(n)}(b)-f^{(n)}(a)\)もほとんどゼロになることが期待されるため、数値計算誤差が限りなく小さくなっていきます。

では変数変換によって端点で減衰する被積分関数に変換することを考えましょう。
高橋-森によって以下の変数変換が提案され、その特性が調べられました[1]。
\(
\displaystyle x=\varphi(t)=\tanh\left(\frac{\pi}{2}\sinh(t)\right)
\)
この変数変換をすることによって端点で”急速に“ゼロに近づく被積分関数に変換されるのです。

この”急速に“の速度はすさまじく、tが大きくなる時、被積分関数は
\(
\displaystyle f(\varphi(t))\varphi'(t)\sim A\exp\{-c\exp(-|t|)\}
\)
の形を持ちます。指数の指数の形で減衰するため、二重指数関数型という名前が付けられているのです。

まぁ、なぜ簡単に計算できるのかと言えば、\(1/\sqrt{x}|_{x\to 0}\)で発散していく時、これよりも早く減衰する関数で割れば収束しますよね、という事です。

※三重、四重指数関数型があるのでは、と考えるのは自然です。しかし、このような関数は(多分ですが)正則であるという条件の下では存在しないことが示されています[3,4]。一重指数関数型数値積分公式は存在します(例えばIMT公式)。

ニュートン・コーツ型数値積分、ガウス求積法は補間型の数値積分公式と呼ばれます。
これらの公式は
与えられたデータ点とデータ点を(重み関数)×(多項式)として補間し、その補間された関数を積分する
というアイデアの元考えられています。
しかし、重み関数は方程式から人間が知っていないとならず、上の形以外の関数を無理やり積分しようとしても積分精度は悪いです。
そのため、問題ごとに合わせたアルゴリズムの変更が必要不可欠になります。

数値計算で気を付けるべき点

端点特異点がある場合の収束

端点特異点がある場合、変数変換後の空間で計算区間が足らなくなり、大きな打切り誤差が発生します。
端点で発散している時に顕著です。
例えば、積分
\(
\displaystyle \int_0^1 \sqrt{x}dx
\)
は高精度(16桁近く一致)の結果を与えますが、
積分
\(
\displaystyle \int_0^1 \frac{1}{\sqrt{x}}dx
\)
は8桁程度の一致しかしません。この原因は打切り誤差です。桁落ちが起きないように処理が必要な問題となります。

積分区間の変更

積分
\(
\displaystyle \int_a^b f(x) dx
\)
を行いたい場合、変数変換
\(
\displaystyle x=\frac{b-a}{2}t+\frac{b+a}{2}
\)
を行うことで区間\([-1,1]\)の積分に置き換えることが出来ます。

広義積分に対する二重指数関数型数値積分

半無限区間\([0,\infty]\)の場合、

\(
\begin{align}
\int_{0}^{\infty} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\exp\left(\frac{\pi}{2}\sinh(t)\right), \\
w_i&=h\frac{\pi}{2}\cosh(ih)\exp\left(\frac{\pi}{2}\sinh(ih)\right)
\end{align}
\)

無限区間\([-\infty,\infty]\)の場合、

\(
\begin{align}
\int_{-\infty}^{\infty} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\sinh\left(\frac{\pi}{2}\sinh(t)\right), \\
w_i&=h\frac{\pi}{2}\cosh(ih)\cosh\left(\frac{\pi}{2}\sinh(ih)\right)
\end{align}
\)

のように変数変換を行うことで二重指数で減衰する関数に変わります。

また、半無限区間の積分の多くの場合は\(exp(-x)\)で減衰します。
これを考慮すると
\(
x=\exp(t-\exp(-t))
\)
という変数変換が有効だということが分かります[1]。本稿では数値計算のアルゴリズムの都合上載せません。

計算の工夫

被積分関数の評価回数を可能な限り減らすため、アルゴリズムを工夫します。

変数変換後の空間で刻み幅\(h\)で求めた数値積分値\(I_h\)は
\(
\displaystyle I_h=h\sum_{i=-N}^N f(\varphi(ih))\varphi'(ih)
\)
と表されます。刻み幅を半分\(h/2\)にすると、分点数は\(4N+1\)になり、
\(
\displaystyle I_{h/2}={h/2}\sum_{i=-2N}^{2N} f(\varphi(ih/2))\varphi'(ih/2)
\)
と表されます。
\(I_h\)を利用して、\(I_{h/2}\)を表現すると、
\(
\displaystyle I_{h/2}=\frac{1}{2}\left\{I_h + h\sum_{i=-N}^{N-1} f(\varphi(ih+h/2))\varphi'(ih+h/2)\right\}
\)
と表されます。

この表式の利点は以下の通りです。
本来、\(I_{h/2}\)を計算するためには\(4N+1\)回右辺を計算しなければなりません。
しかし、\(I_{h}\)と\(I_{h/2}\)の分点の位置が重なっているときがちょうど2N+1点有ります。よって重なっていない残りの\(2N\)点だけを計算してやれば\(I_{h/2}\)を計算できるのです。評価回数が\(4N+1\)回から\(2N\)回に減ったことで単純に計算時間が半分になります。

[adsense2]

fortranプログラム

関数fの数値積分を実行します。
対応しているのは、

・1,2,3次元
・実関数/複素関数(実/複素引数)
・有限区間/半無限区間/無限区間

と

・極座標、球面座標での全空間積分

となっています。
数値計算精度は有効桁数の問題から、8桁程度と思うのが良いと思います。

複素関数の積分は、複素数点と複素数点を直線で結んだ線積分です。
また、半無限区間、無限区間の複素関数の積分では始点\(ca,wa,va\)と角度\(thx,thy,thz\)によって経路を指定することが出来ます。

厳密に全てのルーチンが正しく動作するかは確認していません。
利用する際は連絡先とサイトについてをご覧下さい。

モジュール↓(2000行近くあります)

具体的なプログラム例

積分
\(
\displaystyle \int_0^{5}\sin(\sqrt(x))
\)
を計算。

厳密値
4.3340264879445362
数値計算結果
4.33402648794427

プログラム

全具体例

実装してある計算を全ルーチンを利用したプログラムを書きます。
詳細はコメントを確認してください。

基本的にaを含む変数は積分の始点、bは積分の終点です。

program main
  use DE_integration
  implicit none
  integer::info
  double precision::xa,xb,ya,yb,za,zb,eps,s,thx,thy,thz
  complex(kind(0d0))::ca,cb,wa,wb,va,vb
  double precision,parameter::pi=dacos(-1d0)
  complex(kind(0d0))::cs
 
  double precision::f1,f2,f3,f4,p2,p3,fh1,fi1,fi2,fp2,fi3,fp3
  complex(kind(0d0))::g1,g2,g3,g4,g5,g6,h1,h2,h3,gh1,hh1,&
       gi1,hi1,gi2,hi2,gp2,gi3,hi3,gp3
  external::f1,f2,f3,f4,g1,g2,g3,g4,g5,g6,p2,p3,h1,h2,h3,&
       fh1,gh1,hh1,fi1,gi1,hi1,fi2,gi2,hi2,fp2,gp2,gi3,hi3,gp3,fi3,fp3
 
 
  ! +------- 1D Integration -------+

  ! real f, real x, 1D, finite-range
  xa=0d0;  xb=1d0;  
  eps=1d-12; s=0d0  
  call dde1d(f1,xa,xb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, finite-range
  xa=0d0;  xb=1d0;  
  eps=1d-12; s=0d0  
  call cde1d(g1,xa,xb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, finite-range
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  eps=1d-12; s=0d0  
  call zde1d(h1,ca,cb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x, 1D, semi-infinite
  xa=2d0;  
  eps=1d-12; s=0d0  
  call dde1d_hinf(fh1,xa,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, semi-infinite
  xa=2d0;
  eps=1d-12; s=0d0
  call cde1d_hinf(gh1,xa,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, semi-infinite
  ca=dcmplx(0d0,0d0);  thx=pi/4d0
  eps=1d-12; s=0d0
  call zde1d_hinf(hh1,ca,thx,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x, 1D, infinite
  xa=2d0;  
  eps=1d-12; s=0d0  
  call dde1d_inf(fi1,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, infinite
  xa=2d0;
  eps=1d-12; s=0d0
  call cde1d_inf(gi1,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, infinite
  ca=dcmplx(0d0,0d0);  thx=pi/4d0
  eps=1d-12; s=0d0
  call zde1d_inf(hi1,ca,thx,eps,cs,info)
  write(6,*)info,cs
 
 
  ! +------- 2D Integration -------+
 
  ! real f, real x,y, 2D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0
  eps=1d-12; s=0d0  
  call dde2d(f2,xa,xb,ya,yb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0
  eps=1d-12; s=0d0  
  call cde2d(g2,xa,xb,ya,yb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y, 2D, finite-range
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  wa=dcmplx(1d0,-1d0);  wb=dcmplx(2d0,0d0)  
  eps=1d-12; s=0d0  
  call zde2d(h2,ca,cb,wa,wb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 2D, infinite-range
  eps=1d-8; s=0d0
  call dde2d_inf(fi2,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D, infinite-range
  eps=1d-12; s=0d0  
  call cde2d_inf(gi2,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y, 2D, infinite-range
  ca=dcmplx(0d0,0d0);  thx=0d0
  wa=dcmplx(0d0,0d0);  thy=0d0
  eps=1d-12; s=0d0
  call zde2d_inf(hi2,ca,thy,wa,thx,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 2D polar, infinite-range
  eps=1d-12; s=0d0
  call dde_polar(fp2,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D polar, infinite-range
  eps=1d-12; s=0d0  
  call cde_polar(gp2,eps,cs,info)
  write(6,*)info,cs
 
 
  ! +------- 3D Integration -------+

  ! real f, real x,y,z, 3D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0; za=2d0; zb=3d0
  eps=1d-12; s=0d0
  call dde3d(f3,xa,xb,ya,yb,za,zb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y,z, 3D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0; za=2d0; zb=3d0
  eps=1d-12; s=0d0  
  call cde3d(g3,xa,xb,ya,yb,za,zb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y,z, 3D, finite-range  
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  wa=dcmplx(1d0,-1d0);  wb=dcmplx(2d0,0d0)  
  va=dcmplx(1d0,1d0);  vb=dcmplx(0d0,1d0)
  eps=1d-12; s=0d0  
  call zde3d(h3,ca,cb,wa,wb,va,vb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y,z, 3D, finite-range
  eps=1d-12; s=0d0
  call dde3d_inf(fi3,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y,z, 3D, finite-range
  eps=1d-12; s=0d0  
  call cde3d_inf(gi3,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y,z, 3D, finite-range
  ca=dcmplx(0d0,1d0);  thx=0d0
  wa=dcmplx(1d0,-1d0); thy=pi/2d0
  va=dcmplx(1d0,1d0);  thz=0d0
  eps=1d-12; s=0d0  
  call zde3d_inf(hi3,ca,thx,wa,thy,va,thz,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 3D spherical, infinite-range
  eps=1d-12; s=0d0
  call dde_polar(fp3,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 3D spherical, infinite-range
  eps=1d-12; s=0d0  
  call cde_polar(gp3,eps,cs,info)
  write(6,*)info,cs
 
  stop
end program main

function f1(x)
  implicit none
  double precision,intent(in)::x
  double precision::f1
 
  f1=sin(sqrt(x))
 
  return
end function f1

function g1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::g1
 
  g1=dcmplx(sin(sqrt(x)),exp(-x))
 
  return
end function g1

function h1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::h1
 
  h1=sin(z)
 
  return
end function h1

function fh1(x)
  implicit none
  double precision,intent(in)::x
  double precision::fh1
 
  fh1=exp(-x)
 
  return
end function fh1

function gh1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::gh1
 
  gh1=dcmplx(sqrt(x)*exp(-x),exp(-x))
 
  return
end function gh1

function hh1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::hh1
  double precision,parameter::pi=dacos(-1d0)
 
  hh1=exp(dcmplx(0d0,1d0)*0.5d0*pi*z**2)
 
  return
end function hh1

function fi1(x)
  implicit none
  double precision,intent(in)::x
  double precision::fi1
 
  fi1=exp(-x*x)
 
  return
end function fi1

function gi1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::gi1
 
  gi1=dcmplx(1d0/(1d0+x*x),exp(-x*x))
 
  return
end function gi1

function hi1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::hi1
  double precision,parameter::pi=dacos(-1d0)
 
  hi1=1d0/(1d0+z*z)
 
  return
end function hi1


function f2(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::f2
 
  f2=sin(sqrt(x))*exp(-y)
 
  return
end function f2

function g2(x,y)
  implicit none
  double precision,intent(in)::x,y
  complex(kind(0d0))::g2
 
  g2=dcmplx(sin(sqrt(x))/sqrt(y),exp(-x)*y)
 
  return
end function g2

function h2(z,w)
  implicit none
  complex(kind(0d0)),intent(in)::z,w
  complex(kind(0d0))::h2
 
  h2=sin(z)*w
 
  return
end function h2

function fi2(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::fi2
 
  fi2=1d0/(1d0+x**4+y**4)
 
  return
end function fi2

function gi2(x,y)
  implicit none
  double precision,intent(in)::x,y
  complex(kind(0d0))::gi2
 
  gi2=dcmplx(exp(-y*y-x*x),exp(-y*y-x*x))
 
  return
end function gi2

function hi2(z,w)
  implicit none
  complex(kind(0d0)),intent(in)::z,w
  complex(kind(0d0))::hi2
 
  hi2=exp(-z*z-w*w)
 
  return
end function hi2

function fp2(r,theta)
  implicit none
  double precision,intent(in)::r,theta
  double precision::fp2
 
  fp2=r*exp(-r)*sin(theta)**2
  fp2=r*fp2 ! Jacobian
 
  return
end function fp2

function gp2(r,theta)
  implicit none
  double precision,intent(in)::r,theta
  complex(kind(0d0))::gp2
 
  gp2=dcmplx(0d0,1d0)*r*exp(-r)*sin(theta)**2
  gp2=r*gp2 ! Jacobian
 
  return
end function gp2

function f3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  double precision::f3
 
  f3=sin(x*z)*exp(-y)
 
  return
end function f3

function g3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  complex(kind(0d0))::g3
 
  g3=dcmplx(sin(sqrt(x))/sqrt(y),exp(-x)*y*sqrt(z))
 
  return
end function g3

function h3(z,w,v)
  implicit none
  complex(kind(0d0)),intent(in)::z,w,v
  complex(kind(0d0))::h3
 
  h3=sin(z)*w*v**2
 
  return
end function h3

function fi3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  double precision::fi3
 
  fi3=exp(-x*x-y*y-z*z)
 
  return
end function fi3

function gi3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  complex(kind(0d0))::gi3
 
  gi3=dcmplx(exp(-y*y-x*x-z*z),2d0*exp(-y*y-x*x-z*z))
 
  return
end function gi3

function hi3(z,w,v)
  implicit none
  complex(kind(0d0)),intent(in)::z,w,v
  complex(kind(0d0))::hi3
 
  hi3=exp(-abs(z)**2-abs(w)**2-abs(v)**2)
 
  return
end function hi3

function fp3(r,theta,phi)
  implicit none
  double precision,intent(in)::r,theta,phi
  double precision::fp3
 
  fp3=r*exp(-r)*sin(theta)**2
  fp3=r*r*sin(theta)*fp3 ! Jacobian
 
  return
end function fp3

function gp3(r,theta,phi)
  implicit none
  double precision,intent(in)::r,theta,phi
  complex(kind(0d0))::gp3
 
  gp3=dcmplx(0d0,1d0)*r*exp(-r)*sin(theta)**2
  gp3=r*r*sin(theta)*gp3 ! Jacobian
 
  return
end function gp3

参考文献

[1] H. Takahasi and M. Masatake, “Double Exponential Formulas for Numerical Integration” (1974)
[2]渡辺二太, 二重指数関数型数値積分公式について(1990)

[3] 森 正武、「数値解析における二重指数関数型変換の最適性」
[4] M. Sugihara, Optimality of the double exponential formula-functional analysis approach,
Numer. Math. 75 (1997) 379-395.

数値積分法の入門とDE公式の説明として、
[5] 戸田 英雄, 小野 令美 数値解析における一つの話題

[6] 森 正武「二重指数関数型変換のすすめ」

[7]森 正武, “FORTRAN 77 数値計算プログラミング”、p.174-176岩波書店、(1986)第一刷

プログラミングと数値計算, 数学

畑政義による写像

コメントする

畑政義による写像によって得られる点の集合は、簡単な式で書かれるにもかかわらず、その綺麗さに惹きつけられます
このページでは、畑政義による写像の定義、fortranコード、図を掲載します。

  1. 序
  2. 畑政義による写像の数式
  3. 代表的なパラメータ
  4. パラメータを予想する(2つの複素定数がゼロの場合)
  5. 考察
  6. 参考文献

序

ある日、こんなツイートを見ました。


とても綺麗で感動しました。また、


というツイートに触発されました。自分でも作ってみよう、と。

畑政義による写像の数式

畑政義によって考えられた写像は、簡単な数式で綺麗な図が得られます。

このページでは畑政義の写像で得られる画像の”綺麗さ”を主とします。

畑政義による写像の数式は論文[1]より、
\(
\begin{align}
F_1(z)&=\alpha z+\beta \bar{z} \\
F_2(z)&=\gamma (z-1)+\delta (\bar{z}-1)+1
\end{align}
\)
です(\(\alpha,\beta,\gamma,\delta\)は複素定数,\(\bar{z}\)は\(z\)の複素共役を意味します)。

ある1点に対して写像を行うと、\(F_1,F_2\)によっての2つ点に写像されます。
ある1点からスタートして、\(n\)回作用させていきます。\(n\)回作用させた結果、最後に得られる点の数は\(2^n\)点になります。

具体的にどういう風に計算をすればいいかといいますと
1.初期値\(z_0\)を用意
2.値\(F_1(z_0),F_2(z_0)\)を計算
3.値\(F_1(F_1(z_0)),F_1(F_2(z_0)),F_2(F_1(z_0)),F_2(F_2(z_0))\)を計算
4…
という具合に計算していくのです。そうして得られた最後の複素数の組を複素平面上に打っていくのです。

fortranコードはこちら。重複して計算するのが嫌なので少し工夫しています。

program main
  implicit none
  integer,parameter::N=12
  integer::i
  complex(kind(0d0))::a,b,c,d,z0,h(1:2**N)
 
  a=dcmplx(0.7d0,0.2d0)
  b=dcmplx(0.0d0,0.0d0)
  c=dcmplx(0d0,0d0)
  d=dcmplx(2d0/3d0,0d0)
  z0=dcmplx(1d0,0d0)

  h=dcmplx(0d0,0d0)

  call hatamap(N,a,b,c,d,z0,h)
  do i=1,2**N
     write(10,'(2f10.6)')h(i)
  enddo
 
  stop
end program main

subroutine hatamap(N,a,b,c,d,z0,h)
  implicit none
  integer,intent(in)::N
  complex(kind(0d0)),intent(in)::a,b,c,d,z0
  complex(kind(0d0)),intent(out)::h(1:2**N)

  integer::i,j,k,l,m
  complex(kind(0d0))::z,F,h0(1:2**N)
  external::F
 
  h=dcmplx(0d0,0d0); h0=dcmplx(0d0,0d0)  
 
  h(1)=z0;  h0(1:2**N)=h(1:2**N)
  do j=1,N
     k=1-2**(N-j)
     l=1-2**(N-j+1)
     do i=1,2**j
        if(mod(i,2).eq.1)l=l+2**(N-j+1)
        k=k+2**(N-j)
        m=mod(i+1,2)+1
        h(k)=F(m,a,b,c,d,h0(l))
     enddo
     h0(1:2**N)=h(1:2**N)
  enddo  
 
  return
end subroutine hatamap

function F(n,a,b,c,d,z)
  implicit none
  integer,intent(in)::n
  complex(kind(0d0)),intent(in)::a,b,c,d,z
  complex(kind(0d0))::F
 
  if(n.eq.1)then
     F=a*z+b*conjg(z)
  elseif(n.eq.2)then
     F=c*(z-1d0)+d*(conjg(z)-1d0)+1d0
  else
     write(6,*)"***error"; stop
  endif
 
  return
end function F

さて、プログラムが正しく動いていることを確かめるため、論文のパラメータで写像を再現してみます。
論文[1]では\((\alpha,\beta,\gamma,\delta)\)のセットとして描いています。
論文[1]ではこんな感じに紹介されています。
hata_paper_c

私のプログラムでは、
(左上) \((0.4614+i0.4614, 0, 0.622-i0.196, 0)\)
(右上) \( (0, 0.3+i0.3, 0, 41/50)\)
(左下) \((0, 0.5+i0.5, 0, -0.5+i0.5)\)
(右下) \((0.4614+i0.4614, 0, 0, 0.2896-i0.585)\)
hatamap_reconst_c
となり、ちゃんと正しく計算されていることが分かります。綺麗ですね。

[adsense1]

代表的なパラメータ

さて綺麗な画像を得るために探さなければならないパラメータは\(\alpha,\beta,\gamma,\delta,z_0\)の全部で5つ。各々複素定数なので、実部と虚部で合計10個のパラメータをいじって綺麗な画像を探さなければなりません。これは多いです。
まずは[2],[3],[4]に掲載されているパラメータで計算を行って見ましょう。
ずらっと並べます。
hatamap2_c

hatamap3_c

hatamap4_c

hatamap5_c

hatamap6_c

hatamap1_c

[adsense2]

パラメータを予想する

さて、綺麗な画像として紹介されているのは大体\(\alpha,\beta,\gamma,\delta\)のうち2つがゼロであり、\(0,\pm0.25,\pm 0.5, \pm0.75, \pm1\)のどれかです。この組み合わせだけでも甚大で、
\( _4C_2\times 9^8=258 280 326\)通りの組み合わせがあります。
2億個の画像を見てこれは綺麗、これは違うを判断することはできません。

注記しておきますが、写像を作る際に\(\alpha,\beta,\gamma,\delta\)の絶対値が1以下である必要はありません。
この条件は恐らく指数関数の発散を最低限防ぐためという予想だと思います。

おおよその傾向として以下の推測が出来ました。
\(\kappa, \lambda\)を\(\alpha,\beta,\gamma,\delta\)のゼロではないどれかだとすると、

  • \({\rm Re}{\kappa}={\rm Re}{\lambda}=0\)の時、つまらない画像になる
  • \({\rm Im}{\kappa}={\rm Im}{\lambda}=0\)の時、つまらない画像になる
  • \(\kappa=0\)または\(\lambda=0\)はつまらない画像になる
  • \(\alpha=\beta=0\)または\(\gamma=\delta=0\)はつまらない画像になる

だと感覚的にわかったので、これらを除外します。
さらに、\(0,\pm 0.5, \pm1\)のみを考えて僕自身が綺麗だと思う画像のみを集めてみました。
都合上、\(h=0.5\)と置いています。
また、範囲を固定するために点の撮りうる最大の値を固定してあります。

hatamap001_c3

hatamap002_c3

hatamap003_c3

hatamap004_c3

hatamap005_c3

hatamap006_c3

hatamap007_c3

考察

簡単に分かったことがあります。図の傾向は0.5の数や1の数に密接に関係しているようです。
ここで考えている4つのパラメータは、\(\kappa, \lambda\)のそれぞれ実部、虚部です。
ここで、この4つのパラメータを\((i,j,k,l)\)の組と考えましょう。括弧内の順番は特に関係なく、4つのパラメータの内どれか、を意味するとします。

例えば、初めの方にある、長方形の頂点に点が置かれたような図が得られるとき、4つのパラメータの組は\((0,0,\pm 1,\pm h)\)の組み合わせが多いことに気が付きます。
h100_c

続いて、これぞフラクタルだ、と思わせる形は\((\pm h, \pm h, \pm h, \pm h)\)、
hhhh_c

最後の方の紋様のような模様では\((\pm 1, \pm 1, \pm 1, \pm 1)\)で埋められているようです。
1111_c

恐らく、フラクタル的な雰囲気を持つ画像は、\(\kappa, \lambda\)の各々の値が有限で、絶対値が1とか、そんな時に出てくるのではないでしょうか。

今回調べた範囲のは、\(\alpha,\beta,\gamma,\delta\)のうち2つがゼロ、しかも0.5刻みしか許さないという非常に強い条件を入れた範囲です。
これだけでも膨大な数になることが、上の画像の量を見ただけで分かるでしょう。

本来は\(\alpha,\beta,\gamma,\delta\)がゼロではなくて良く、0.5刻みである必要もなく、実数ですから、今回調べたのはごくごく一部です。

綺麗な画像の定義があいまいなのもいけませんね。どうにかしてみたいものです。

また、同じく@snowy_treeさんが、


という綺麗なパラメータを集めたgifを公開してくれました。凄く綺麗なパラメータがたくさんあります！
とんでもなく綺麗な写像を見てみたいものです。

参考文献

[1]Dimensión, Análisis Multifractal y Aplicaciones
内の「3. Fractals in Mathematics (Hata)」
※中を見ると掲載されている元の論文はPatterns and Waves pp.259-278(1986)という論文のようです。
畑さんは、写像を拡張し統一的に扱おうとしたようです。なので、畑さんがこの写像を1から全て作った、というのは語弊があるかと思います。

[2]カオス #3,#6,#10,#11,#12,#13,#14,#15,#27畑政義写像 -主にコーディング

Web上で畑政義による写像が書けるサイトがありますので紹介いたします。
[2]畑政義写像で遊ぼう -救済に紹介されています。
(Play with Hata-map / 畑政義写像で遊ぼう)

[3]畑政義写像(1) -閃光的網站・弛緩複合体 Blog
※1. ただし、若干異なっています。
n回繰り返した時のみの写像の集まりのはずですが、どうやらn回に至るまでの点もプロットしているようです。n+1回繰り返しても似た画像しか出てこないのでまぁ、問題ないと思いますが…。

※2. また、どちらも上下が逆さまになっています。時々この上の二つで作成した値を私のプログラムに入れても再現できない場合があります。これに起因する問題なのか、別の要因なのかははっきりしません。

プログラミングと数値計算, 数学

平均、分散、標準偏差

コメントする

(i) 全データが存在する場合

　例１ ： ４０人クラスにおいてテスト点数の解析を行いたい(４０人全員の点数を知っている)。
　例２ ： 過去１年間、365日で体重がどのくらい増減したかの解析をしたい(365日分のデータが存在する)。
平均\(m\)
\(
\displaystyle m=\frac{1}{N}\sum_{i=1}^N x_i
\)
分散\(v\)
\(
\displaystyle v=\frac{1}{N}\sum_{i=1}^N (x_i-m)^2
\)
標準偏差\(\sigma\)
\(
\displaystyle \sigma=\sqrt{v}
\)

平均値から\(\pm 66\%\)の誤差に収まっている範囲は、
平均m ± 標準偏差σ　 と書ける。
例1 のクラスの点数の場合は、
平均m ± 標準偏差σ [点]　　となる。

(ii)一部分のデータしか存在しない場合

例１ ： ４０人クラスにおいて２３人しかテストの点数を知らないが、その２３人からクラス全体の点数の解析を行いたい（２３人の点数は分かる。残り１７人の点数は不明）。
　例２ ： 過去１年間、365日で体重の増減を解析したいが、200日分のデータしかない。だけど１年の解析を行いたい（165日分のデータは不明、紛失した）。
　　　　例３ ： 物理量の測定（測定回数が1000回だとしても、物理量は無限回の測定をしなければ真値は不明だ、と考える。）

この場合、平均、分散、標準偏差は
平均\(m’\)
\(
\displaystyle m’=\frac{1}{N}\sum_{i=1}^N x_i
\)
分散\(v’\)
\(
\displaystyle v’=\frac{1}{N}\sum_{i=1}^N (x_i-m’)^2
\)
標準偏差\(\sigma’\)
\(
\displaystyle \sigma’=\sqrt{v’}
\)
平均値から\(\pm 66\%\)の不確かさに収まっている範囲は、
平均m’ ± 標準偏差σ’　 と書ける。
例1 のクラスの点数の場合は、
平均m’ ± 標準偏差σ’ [点]　　となる。

誤差と不確かさ　は違うので注意しましょう。
誤差 → 完璧なデータのときに使う。
不確かさ → 本当の値を知らないときに使う。

プログラミングと数値計算, 数学

ディガンマ関数の数値計算

4件のコメント

ディガンマ関数\(\psi(z)\)は

\(
\displaystyle \psi(z)=\frac{d}{dz} \ln\Gamma(z)=\frac{\Gamma'(z)}{\Gamma(z)}
\)
として定義されます。

数値計算でこれを求めます。
用いる式は、
ディガンマ関数の漸近展開式[1],
\(
\displaystyle \psi(z)\sim \ln z -\frac{1}{2z}-\sum_{n=1}^{\infty} \frac{B_{2n}}{2nz^{2n}},\;\; z\to\infty (|{\rm arg} z| \lt \pi )
\)
ここで\(B_{2n}\)はベルヌーイ数[2]
と
漸化式[1]
\(
\displaystyle \psi(z+1)=\psi(z)+\frac{1}{z}
\)
相反公式[1]
\(
\psi(1-z)=\psi(z)+\pi \cot(\pi z)
\)
を利用します。

具体的には領域を4つに分けて、

  1. \(10 \lt {\rm Re}(z)\),      漸近展開
  2. \(1 \lt {\rm Re}(z) \lt 10\),    漸近展開+漸化式
  3. \(-9 \lt {\rm Re}(z) \lt 1\),    漸近展開+漸化式+相反公式
  4. \({\rm Re}(z) \lt -9\),      漸近展開+相反公式

を利用しています。

program main
  implicit none
  integer::i,j
  double precision::x,y
  complex(kind(0d0))::z,digamma
  external::digamma
 
  do i=0,200
     x=-10.05d0+i*0.1d0
     do j=0,200
        y=-10.05d0+j*0.1d0
        z=dcmplx(x,y)
        write(10,'(4e20.10e2)')x,y,dble(digamma(z)),dimag(digamma(z))
     enddo
     write(10,*)
  enddo
 
  stop
end program main

!--------------------------

function digamma(z)
  ! digamma function for complex plane
  ! sikinote, http://slpr.sakura.ne.jp/qp/
  ! author : sikino
  ! date : 2016/07/04 (yyyy/mm/dd)
  ! date : 2016/07/07 (yyyy/mm/dd)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::digamma
 
  integer::i,j,n,m,as,key
  integer,parameter::nmax=100
  double precision::eps=1d-13
  double precision::pi=4d0*atan(1d0)
  complex(kind(0d0))::c,w,t,ds,s
  double precision::B
  external::B

  as=10

  if(dble(z).ge.dble(as))then
     w=z
     key=0
  elseif(dble(z).ge.1d0)then
     w=z+dble(as)
     key=1
  elseif(dble(z).ge.-dble(as-1))then
     w=1d0-z+dble(as)
     key=2
  else
     w=1d0-z
     key=3
  endif
 
  ! Asymptotic expansion
  s=log(w)-0.5d0/w
  do n=1,nmax
     m=2*n
     ds=B(m)/(m*w**m)
     s=s-ds
     if(abs(ds)/abs(s).le.eps)exit
  enddo
  if(n.eq.nmax+1)then
     write(6,*)"Worning! did not converge at digamma"
  endif
 
  if(key.eq.1.or.key.eq.2)then
     !Recurrence relation
     do i=1,as
        s=s+1d0/(1d0-w)
        w=w-1d0
     enddo
  endif
  if(key.eq.2.or.key.eq.3)then
     ! Reflection Formula
     s=s-pi/tan(pi*z)
  endif

  digamma=s
 
  return
end function digamma
!--------------------------
function B(n)
  ! 1st bernoulli number B(n), 0 <= n <= 100
  ! from wolfram alpha
  !    " Table[bernoulli b(i), {i,0,100}] with 17 digits "
  implicit none
  integer,intent(in)::n
  double precision::B,A(0:100)
 
  if(n.lt.0.or.n.gt.100)then
     write(6,*)"Not expected n at bernoulli number function. "
     write(6,*)"program stop"
     stop
  endif
 
  A(0:100)=(/1d0 &
       ,  -0.50000000000000000d0, 0.16666666666666667d0, 0d0, -0.033333333333333333d0, 0d0 &
       ,  0.023809523809523810d0, 0d0, -0.033333333333333333d0, 0d0,  0.075757575757575758d0, 0d0 &
       ,  -0.25311355311355311d0, 0d0,    1.1666666666666667d0, 0d0,   -7.0921568627450980d0, 0d0 &
       ,    54.971177944862155d0, 0d0,   -529.12424242424242d0, 0d0,    6192.1231884057971d0, 0d0 &
       ,   -86580.253113553114d0, 0d0,   1.4255171666666667d+6, 0d0,  -2.7298231067816092d+7, 0d0 &
       ,   6.0158087390064237d+8, 0d0, -1.5116315767092157d+10, 0d0,  4.2961464306116667d+11, 0d0 &
       , -1.3711655205088333d+13, 0d0,  4.8833231897359317d+14, 0d0, -1.9296579341940068d+16, 0d0 &
       ,  8.4169304757368262d+17, 0d0, -4.0338071854059455d+19, 0d0,  2.1150748638081992d+21, 0d0 &
       , -1.2086626522296526d+23, 0d0,  7.5008667460769644d+24, 0d0, -5.0387781014810689d+26, 0d0 &
       ,  3.6528776484818123d+28, 0d0, -2.8498769302450882d+30, 0d0,  2.3865427499683628d+32, 0d0 &
       , -2.1399949257225334d+34, 0d0,  2.0500975723478098d+36, 0d0, -2.0938005911346378d+38, 0d0 &
       ,  2.2752696488463516d+40, 0d0, -2.6257710286239576d+42, 0d0,  3.2125082102718033d+44, 0d0 &
       , -4.1598278166794711d+46, 0d0,  5.6920695482035280d+48, 0d0, -8.2183629419784576d+50, 0d0 &
       ,  1.2502904327166993d+53, 0d0, -2.0015583233248370d+55, 0d0,  3.3674982915364374d+57, 0d0 &
       , -5.9470970503135448d+59, 0d0,  1.1011910323627978d+62, 0d0, -2.1355259545253501d+64, 0d0 &
       ,  4.3328896986641192d+66, 0d0, -9.1885528241669328d+68, 0d0,  2.0346896776329074d+71, 0d0 &
       , -4.7003833958035731d+73, 0d0,  1.1318043445484249d+76, 0d0, -2.8382249570693707d+78/)
 
  B=A(n)
 
  return
end function B

[adsense2]

参考文献

[1]Abramowitz, M. and I.A. Stegun, {\it Handbook of Mathematical Functions}, Dover Publications, Ninth Printing, 1970., P.258-259, Eq.6.3.18(漸近展開式), Eq.6.3.5(漸化式), Eq.6.3.7(相反公式)
[2]Bernoulli Number -wolfram mathworld

プログラミングと数値計算, 数学

ベルヌーイ数の数値計算

コメントする

これがベルヌーイ数を求める最速のアルゴリズムだと思います。
直接代入。これに限ります。

n=0~100までサポートしています。
実用的にはこの範囲で足りるんじゃないかなと思います。

値はwolfram alphaより。
Table[bernoulli b(i), {i,0,100}] with 17 digits

program main
  implicit none
  integer::i
  double precision::B
  external::B
 
  do i=0,30
     write(6,'(i5,f30.17)')i,B(i)
  enddo
 
  stop
end program main
   
function B(n)
  ! 1st bernoulli number B(n), 0 <= n <= 100
  ! from wolfram alpha
  !    " Table[bernoulli b(i), {i,0,100}] with 17 digits "
  implicit none
  integer,intent(in)::n
  double precision::B
  double precision::A(0:100)
 
  if(n.lt.0.or.n.gt.100)then
     write(6,*)"Not expected n at bernoulli number function. "
     write(6,*)"program stop"
     stop
  endif
   
  A(0:100)=(/1d0 &
       ,  -0.50000000000000000d0, 0.16666666666666667d0, 0d0, -0.033333333333333333d0, 0d0 &
       ,  0.023809523809523810d0, 0d0, -0.033333333333333333d0, 0d0,  0.075757575757575758d0, 0d0 &
       ,  -0.25311355311355311d0, 0d0,    1.1666666666666667d0, 0d0,   -7.0921568627450980d0, 0d0 &
       ,    54.971177944862155d0, 0d0,   -529.12424242424242d0, 0d0,    6192.1231884057971d0, 0d0 &
       ,   -86580.253113553114d0, 0d0,   1.4255171666666667d+6, 0d0,  -2.7298231067816092d+7, 0d0 &
       ,   6.0158087390064237d+8, 0d0, -1.5116315767092157d+10, 0d0,  4.2961464306116667d+11, 0d0 &
       , -1.3711655205088333d+13, 0d0,  4.8833231897359317d+14, 0d0, -1.9296579341940068d+16, 0d0 &
       ,  8.4169304757368262d+17, 0d0, -4.0338071854059455d+19, 0d0,  2.1150748638081992d+21, 0d0 &
       , -1.2086626522296526d+23, 0d0,  7.5008667460769644d+24, 0d0, -5.0387781014810689d+26, 0d0 &
       ,  3.6528776484818123d+28, 0d0, -2.8498769302450882d+30, 0d0,  2.3865427499683628d+32, 0d0 &
       , -2.1399949257225334d+34, 0d0,  2.0500975723478098d+36, 0d0, -2.0938005911346378d+38, 0d0 &
       ,  2.2752696488463516d+40, 0d0, -2.6257710286239576d+42, 0d0,  3.2125082102718033d+44, 0d0 &
       , -4.1598278166794711d+46, 0d0,  5.6920695482035280d+48, 0d0, -8.2183629419784576d+50, 0d0 &
       ,  1.2502904327166993d+53, 0d0, -2.0015583233248370d+55, 0d0,  3.3674982915364374d+57, 0d0 &
       , -5.9470970503135448d+59, 0d0,  1.1011910323627978d+62, 0d0, -2.1355259545253501d+64, 0d0 &
       ,  4.3328896986641192d+66, 0d0, -9.1885528241669328d+68, 0d0,  2.0346896776329074d+71, 0d0 &
       , -4.7003833958035731d+73, 0d0,  1.1318043445484249d+76, 0d0, -2.8382249570693707d+78/)
 
  B=A(n)
 
  return
end function B

B(256)まで増やしたのはこちら。

一応逐次的に求める事が出来るアルゴリズムがありまして実装しましたが、全然良くないです[1]。
4倍精度演算を使って\(B_{30}\)を求めようとするとたったの8桁しか一致しません。
倍精度演算のみだと\(B_{14}\)を求めようとするとたったの6桁しか一致しません。

参考文献

[1]Algorithmic description -Bernoulli number, wikipedia

[adsense2]

プログラミングと数値計算, 数学

リーマンのゼータ関数の数値計算コード(複素平面)

8件のコメント

ここではゼータ関数を数値計算で得るためのアルゴリズムとfortran90によるプログラムを載せています。

目次

  1. ゼータ関数の級数表示
  2. プログラム(fortran90)
  3. 計算データ
  4. グラフ
  5. 非自明、自明なゼロ点
  6. ゼータ関数の導関数とプログラム

ゼータ関数の計算方法

ゼータ関数は通常
・\(\rm{Re}(s)\gt 1\)に対しては直接計算を行い、\(\rm{Re}(s)\le 1\)の領域に対しては相反公式
\(
\displaystyle \zeta(s)=2^{s}\pi^{s-1}\sin\left(\frac{\pi s}{2}\right)\Gamma(1-s)\zeta(1-s)
\)
を使って求める事が推奨されます。
ガンマ関数を使うことを避けたい場合、全複素平面(\(s\ne 1\))で有効なゼータ関数の表示
\(
\displaystyle
\zeta(s)=\frac{1}{1-2^{1-s}}
\sum_{n=0}^{\infty} \frac{1}{2^{n+1}}
\sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
\)
を使うのが良いでしょう。

ゼータ関数の級数表示

ゼータ関数は
\(
\displaystyle
\zeta(s)=\sum_{n=1}^{\infty} \frac{1}{n^s} ,~~~(1\lt s)
\)
と定義されます。が、これはlogの収束であるため非常に収束が遅いです。
Euler-Knopp変換
\(\displaystyle
\sum_{n=0}^\infty a_n = \sum_{n=0}^{\infty} \frac{1}{(1+q)^{n+1}}
\sum_{j=0}^{n}{n \choose j} q^{n-j}a_j
\)
(q=1のときEuler変換)によって早く収束する級数に変換されます[4]。
より詳しい証明は[5]。

Riemann Zeta Function — from Wolfram MathWorld, Eq.(21)
によると、全領域(s=1以外)で収束するリーマンのゼータ関数の級数表示は、
\(
\displaystyle
\zeta(s)=\frac{1}{1-2^{1-s}}
\sum_{n=0}^{\infty} \frac{1}{2^{n+1}}
\sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
\)
と書けます。もう一つ形式があるのですが、そちらは収束が遅そうな形式だったので用いていません。

数値計算的には2の乗数とかはあまり計算したくありません。そこで、漸化式を考えます。
\(
\begin{align}
\zeta(s) &=\sum_{n=0}^{\infty} b_n \\
b_n &= \sum_{k=0}^{n}a_k^n \\
a_k^n &= \frac{1}{1-2^{1-s}}\frac{1}{2^{n+1}}{n \choose k}\frac{(-1)^k}{(k+1)^s}
\end{align}
\)
と\(b_n, a_k^n\)を定義します。\(a_k^n\;\; (n\ge k)\)を計算するために漸化式を考えると

\(
\begin{align}
a_k^n&=\frac{1}{2}\left(\frac{n}{n-k}\right)a_{k}^{n-1} \\
a_k^n&=\frac{k-n-1}{k}\left(\frac{k}{k+1}\right)^s a_{k-1}^n \\
a_0^0&=\frac{1}{2}\frac{1}{1-2^{1-s}}
\end{align}
\)

と導く事が出来ます。
計算手順としては、

zeta_alg_c

としていけばよいでしょう。
計算量を考えて赤の手順を少なくするようにしています。

fortran90によるプログラム

本当のゼータ関数の値との相対誤差\(\delta\)は約
\(\delta\lt 3\times 10^{-14} (Re(s)\lt 1)\),
\(\delta\lt 3\times 10^{-15} (1 \lt Re(s)) \)
負の領域の方が精度が悪いのは、相反公式から来ています。
相反公式のみを4倍精度演算で行えば
\(\delta\lt 7\times 10^{-15} (Re(s) \lt 1) \)
にまで精度が上がります。目立って精度はあがりません。

2017/01/12に大幅な更新をしました。
以前のものと比べ変更点は、
・4倍精度を使用していない
・計算速度は以前の54倍
になっているので、以前のものを使っている方は更新することを推奨します。
ガンマ関数に関しては複素ガンマ関数の数値計算をご参照ください。

function zeta(s)
  implicit none
  complex(kind(0d0)),intent(in)::s
  complex(kind(0d0))::zeta
  !
  !Checked : ifort   (success : version 16.0.2 20160204)
  !        : gfortran(success : version 4.8.4)
  !        : gfortran(success : version 4.4.7)
  !
  !http://slpr.sakura.ne.jp/qp/
  ! Author : sikino
  ! date: 2016/06/02, (yyyy/mm/dd)
  !     : 2016/06/03
  !     : 2016/06/04
  !     : 2016/06/05
  !     : 2016/12/25
  !     : 2017/01/12
  !     : 2017/01/14
  !
  integer::n,k
  integer,parameter::nmax=1000
  double precision::nq,nq1,ln1,ln2
  double precision,parameter::eps=1d-15
  double precision,parameter::pi=3.14159265358979324d0
  complex(kind(0d0))::ds,is,d,s0,tc,tt,qzeta
  complex(kind(0d0))::a(0:Nmax),b(0:nmax),cgamma
  external::cgamma
  !complex*32::zgamma
  !external::zgamma

  ds=dble(s); is=imag(s)
  if(ds.eq.1d0.and.is.eq.0d0)then
     zeta=1d+300
  elseif(ds.eq.0d0.and.is.eq.0d0)then
     zeta=-0.5d0
  else
     if(real(s).lt.1d0)then
        s0=1d0-s
     else
        s0=s
     endif
     
     tt=1d0/(1d0-2d0**(1d0-s0))    
     qzeta=0.5d0*tt

     a=cmplx(0d0,0d0,kind=8); b=cmplx(0d0,0d0,kind=8)
     b(0)=0.25d0*tt
     do n=1,nmax
        nq=n
        nq1=nq+1d0

        a(0)=b(0)
        do k=1,n-1
           a(k)=b(k)*nq/(n-k)
        enddo
        !----
        ln2=log(dble(n)/dble(n+1))        
        if(n.eq.1)then
           d=s0*ln2
        else
           d=d*ln2/ln1
        endif
        a(n)=-a(n-1)*exp(d)/nq
        ln1=ln2
        !----

        tc=cmplx(0d0,0d0,kind=8)
        do k=0,n
           tc=tc+a(k)
        enddo

        qzeta=qzeta+tc

        if(abs(qzeta).ge.1d0)then
           if(abs(tc/qzeta).le.eps)exit
        else
           if(abs(tc).le.eps)exit
        endif

        b=0.5d0*a
     enddo
     if(n.eq.nmax+1)write(6,*)" ***worning, didn't converge at zeta"

     if(real(s).lt.1d0)then
        zeta=(2d0**s)*pi**(-s0)*sin(0.5d0*pi*s)*qzeta*cgamma(s0)
        !zeta=cmplx((2q0**s)*acos(-1q0)**(-s0)*sin(0.5q0*acos(-1q0)*s)*qzeta*zgamma(s0*1q0),kind=8)
    else
        zeta=qzeta
     endif
  end if

  return
end function zeta

function cgamma(z)
  ! sikinote (http://slpr.sakura.ne.jp/qp)
  ! author : sikino
  ! date   : 2017/01/09
  !          2017/01/10
 
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::cgamma
 
  integer::i,n,M
  double precision,parameter::pi=3.14159265358979324d0
  double precision,parameter::e1=0.3678794411714423216d0
  double precision,parameter::ln2pi2=0.91893853320467274d0
  double precision,parameter::sq2pi=2.5066282746310005d0
  double precision::a(1:30),AB(1:14),dz,iz
  complex(kind(0d0))::q,s,w,r,q1,q2

  dz=dble(z)
  iz=imag(z)

  if(dz.eq.nint(dz).and.iz.eq.0d0)then
     if(dz.gt.0d0)then
        s=dcmplx(1d0,0d0)
        n=nint(dz)
        do i=2,n
           s=s*i
        enddo
     else
        s=1d+300
     endif
     cgamma=s
  else
     q=z
     if(dz.lt.0d0)q=1d0-z

     if(abs(iz).lt.1.45d0)then
        ! For x=arb, |y|\lt 1.45
        n=int(dble(q))
        s=1d0
        do i=1,n
           s=s*(q-i)
        enddo
        w=q-dble(n)

        a(1) = 1d0
        a(2) = 0.57721566490153286d0
        a(3) =-0.65587807152025388d0
        a(4) =-0.42002635034095236d-1
        a(5) = 0.16653861138229149d0
        a(6) =-0.42197734555544337d-1
        a(7) =-0.96219715278769736d-2
        a(8) = 0.72189432466630995d-2
        a(9) =-0.11651675918590651d-2
        a(10)=-0.21524167411495097d-3
        a(11)= 0.12805028238811619d-3
        a(12)=-0.20134854780788239d-4
        a(13)=-0.12504934821426707d-5
        a(14)= 0.11330272319816959d-5
        a(15)=-0.20563384169776071d-6
        a(16)= 0.61160951044814158d-8
        a(17)= 0.50020076444692229d-8
        a(18)=-0.11812745704870201d-8
        a(19)= 0.10434267116911005d-9
        a(20)= 0.77822634399050713d-11
        a(21)=-0.36968056186422057d-11
        a(22)= 0.51003702874544760d-12
        a(23)=-0.20583260535665068d-13
        a(24)=-0.53481225394230180d-14
        a(25)= 0.12267786282382608d-14
        a(26)=-0.11812593016974588d-15
        a(27)= 0.11866922547516003d-17
        a(28)= 0.14123806553180318d-17
        a(29)=-0.22987456844353702d-18
        a(30)= 0.17144063219273374d-19

        M=int(11.3*abs(w)+13)
        if(M.gt.30)M=30
        r=a(M)
        do i=M-1,1,-1
           r=r*w+a(i)
        enddo
       
        cgamma=s/(r*w)
     else
        ! For x=arb, |y|\gt 1.45
        s=1d0
        if(abs(q).lt.9d0)then
           do i=0,7
              s=s*(q+i)
           enddo
           s=1d0/s
           q=q+8d0
        endif

        AB(1) = 0.83333333333333333d-1
        AB(2) =-0.27777777777777778d-2
        AB(3) = 0.79365079365079365d-3
        AB(4) =-0.59523809523809524d-3
        AB(5) = 0.84175084175084175d-3
        AB(6) =-0.19175269175269175d-2
        AB(7) = 0.64102564102564103d-2
        AB(8) =-0.29550653594771242d-1

        q1=1d0/q
        q2=q1*q1
       
        r=AB(8)
        do i=7,1,-1
           r=r*q2+AB(i)
        enddo
        cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1)
     endif
     if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z))
  endif

  return
end function cgamma

ifort version 16.0.2 20160204, 12.1.0 20110811
gfortran version 4.8.4
で動くことを確かめています。

こちらに4倍精度用のゼータ関数のコードを置いておきます。
精度は何点かで試しましたが、30~32桁一致するくらいです。

ゼータ関数のデータ

計算結果のデータも置いておきます。
(11MB)riemanndata.tar.gz
(9MB)riemanndata.zip
それぞれ、中に
“riemannzeta_high_resolution.txt”
“riemannzeta_low_resolution.txt”
が入っています。
highの方は1列目から4列目にかけて
\({\rm Re}(s):[-20\sim 20],{\rm Im}(s):[-50\sim 50],{\rm Re}(\zeta(s)),{\rm Im}(\zeta(s))\)、0.1ずつ、小数点以下13桁で出力したもので、
lowの方は1列目から4列目にかけて
\({\rm Re}(s):[-20\sim 20],{\rm Im}(s):[-50\sim 50],{\rm Re}(\zeta(s)),{\rm Im}(\zeta(s))\)、0.2ずつ、小数点以下5桁で出力したものです。

x,y,z軸を\({\rm Re}(s),{\rm Im}(s),{\rm Re}(\zeta(s))\)にとりグラフにすると、
riemann_zeta_c
というグラフが得られます。

\({\rm Re}(s)\)が負で\({\rm Im}(s)\)が小さいところで発散が起こるようですので、注意をして使ってください。
何故倍精度でおかしい問題が起こっていたのか？今回の場合はオーバーフローでもなく、アンダーフローでもありません。
原因は桁落ちです。和と差によって激しく桁落ちしてしまうのです。

また、倍精度よりも計算が早いです。なぜならば目標精度までスムーズに収束してくれるからです。

相反公式を用いたので上記問題は発生しません。

[adsense1]

ゼータ関数の実部、虚部、絶対値のグラフ

ゼータ関数の実部、虚部、絶対値のグラフは以下の通りになります。

実部\({\rm Re}(\zeta(s))\)

zeta_real_quad
カラーマップ(横軸:\({\rm Re}(s)\),縦軸\(:{\rm Im}(s)\),カラー:\({\rm Re}(\zeta(s))\))
map_riemann_real

虚数\({\rm Im}(\zeta(s))\)

zeta_imag_quad_c
カラーマップ(横軸:\({\rm Re}(s)\),縦軸\(:{\rm Im}(s)\),カラー:\({\rm Im}(\zeta(s))\))
map_riemann_imag

絶対値\(|\zeta(s)|\)

zeta_abs_quad_c
カラーマップ(横軸:\({\rm Im}(s)\),縦軸\(:{\rm Re}(s)\),カラー:\(|\zeta(s)|\))

対数プロット(横軸:\({\rm Re}(s)\),縦軸\(:{\rm Im}(s)\),カラー:\((\log{10}(|\zeta(s)|\))
map_riemann_logabs

（おまけ）ゼータ関数のニュートン写像(詳しくはリーマンのゼータ関数の複素力学系へ)

非自明なゼロ点

上記の計算方法(4倍精度のプログラム)で、リーマンの非自明なゼロ点を考えます。
「非自明なゼロ点」はゼータ関数の実数軸上には無いゼロ点のことです。
リーマン予想とは
ゼータ関数\(\zeta(s)\)の「非自明なゼロ点」が\(s=\frac{1}{2}+ix\)(\(x\)は実数)上に全てある
というものです。\(\zeta(s)=0\)が満たされるときはゼータ関数の実部も虚部もゼロになるので、簡単に言えばゼータ関数の絶対値がゼロになる\(|\zeta(s)|=0\)と言うことです。
実際に\(|\zeta(s)|\)を可視化してみますとこうなります。
quadzeta_c

自明なゼロ点

僕は他のページでは自明な点に関する図を見たことがありません。非常に小さいからでしょう。
小さいものを見るには対数スケールで見るのが良いです。実際に書かせてみます。
quadzeta_log2_c
確かに小さいですね。対数にしないと明らかに見えません。

絶対値\(|\zeta(s)|\)を上の図のように表示させるgnuplotスクリプトはこちら。

set xr[-15:10]
set yr[-40:40]
set zr[0:2]
set view 69, 73, 2, 0.3
set view  equal xy
set ticslevel 0
   
set palette defined (-9 "purple",-6 "blue",-3 "turquoise", 0 "#f5f5f5", 3 "gold", 6 "orange", 9 "red")
set colorbox horiz user origin .1,.1 size .8,.04
set pm3d at b

splot "fort.11" u 1:2:(sqrt($3**2+$4**2)) every :3 w l lw 0.1

[adsense2]

クリティカルライン上での計算

クリティカルラインとは
\(
\displaystyle s=\frac{1}{2}+ix
\)
の線です。リーマン予想は、非自明なゼロ点が全てこの上にあるか？という問いかけであり、この計算と繋がっています。

倍精度型のアルゴリズム(zeta)では、\(x\sim 400\)程度まで計算することが出来ます。
これは、

相反公式の
\(
\displaystyle \sin(\frac{\pi}{2}s)
\)
の値がオーバーフローを起こすためです。

倍精度演算では約\(10^{306}\)まで扱えるので、
\(
s=\sin^{-1}(10^{306})\frac{2}{\pi}\sim 1+i449
\)
4倍精度では約\(10^{4932}\)まで扱えるので、
\(
s=\sin^{-1}(10^{4932})\frac{2}{\pi}\sim 1+i7230
\)
となります。
しかし、実際に計算してみると倍精度では確かに\(x=450\)程度まで計算可能ですが、4倍精度では\(x=1600\)程度までしか正しい値を返してくれません。どこかで何かが起こっているはずですが、解析はしてません。

実際に試してみますと、
倍精度では

となります。

ゼータ関数の導関数とプログラム

ゼータ関数の導関数も求めましょう。
導関数はゼータ関数の式をそのまま微分すればいいです。なぜなら、正則関数として解析接続されているからです。
実際に微分しますと、
\(
\displaystyle
\zeta'(s)=\frac{1}{(2^s-2)^2}\sum_{n=0}^{\infty} \frac{1}{2^{n+1}} \sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
\left[ -2^s \{(2^s-2)\ln(k+1)+\ln4\} \right]
\)
変形して以下のように
\(
\begin{align}
\zeta'(s) &=\sum_{n=0}^{\infty} b_n \\
b_n &= \sum_{k=0}^{n}a_k^n \\
a_k^n &= (-1)^k {n \choose k} \left[ \frac{-2^s (k+1)^{-s}\{(2^s-2)\ln(k+1)+\ln4\}}{(2^s-2)^2} \right]
\end{align}
\)
となります。同じように漸化式を求めてやれば、
\(
\begin{align}
a_k^n&=\frac{1}{2}\left(\frac{n}{n-k}\right)a_{k}^{n-1} \\
a_k^n&=-\frac{k-n-1}{k}\left(\frac{k}{k+1}\right)^s\frac{(2^s-2)\ln(k+1)+\ln 4}{(2^s-2)\ln(k)+\ln 4} a_{k-1}^n \\
a_0^0&=\frac{1}{2}\frac{-2^s \ln 4}{(2^s-2)^2}
\end{align}
\)

として求められます。
実際の数値計算では、2番目の漸化式は余り有能ではないかもしれません。
ゼータ関数を求めるときと同じようにやるならば、必要なのは\(a_n^n\)を\(a_{n-1}^n\)から求める作業です。漸化式を用いても直接求めても計算量はさほど変わらない気がします。

導関数の相反公式

相反公式ももちろん存在します[3]。
解析接続されているので、元の相反公式の左辺、右辺をそのまま微分することで求められます。
\(
\begin{align}
-\zeta'(1-s)&=2(2\pi)^{-s}\left\{x\cos\left(\frac{\pi s}{2}\right)+y\sin\left(\frac{\pi s}{2}\right)\right\}\Gamma(s)\zeta(s) \\
&\hspace{3em} +2(2\pi)^{-s}\cos\left(\frac{\pi s}{2}\right)\left\{\Gamma(s)\zeta(s)\right\}’ \\
& x = -\ln(2\pi), \;\; y=-\frac{\pi}{2}
\end{align}
\)

\(
\Gamma'(s)=\Gamma(s)\psi(s)
\)
ここで\(\psi(z)\)はディガンマ関数

確かめはwolfram alphaを用いて以下のコマンドで行いました。
ReplaceAll[D[Zeta[x], x], {x -> -6+2i}]

fortran90によるプログラム

fortran90によるプログラムはこんな感じになるかと思います。
相反公式で、ゼータ関数、倍精度ガンマ関数とディガンマ関数を用います。
ゼータ関数は上の、倍精度ガンマ関数とディガンマ関数はつけておきました。詳しくはディガンマ関数の数値計算をご覧ください。
ここに記述しているのは倍精度のみです。

注釈

他のページでの例えば
リーマンのゼータ関数で遊び倒そう (Ruby編)
等では、漸化式を使わないで直接計算しています。このため、余り項の数\(n\)の数を増やす事が出来ないこと、倍精度と思われるため、桁落ちの問題が回避できない問題が起こります。
実際、私も直接計算するコードをfortran90で書きましたが、項の数を増やせず収束しない事、また計算時間が2~3倍余計にかかります。

参考文献
Riemann Zeta Function — from Wolfram MathWorld, Eq.(21)
draw a horizontal gradient somewhere at the bottom of the graph. -ColorBox

[3]Tom M. Apóstol “Formulas for Higher Derivatives of the Riemann Zeta Function”, Math. of. Comp vol 44, 169, 1985, pp. 223-232.

[4]Madhava-Gregory-Leibniz級数の収束の加速(収束の加速1)

[5]リーマンゼータ関数の級数表示による解析接続 インテジャーズ

プログラミングと数値計算, 数学

特殊関数

3件のコメント

以下に示すのは

・ラゲール陪多項式(Associated Laguerre Polynomials):\(L_n^{(\alpha)}(x)\)
・ルジャンドル陪多項式(Associated Legendre Polynomials):\(P_l^m(x)\)
・ヤコビ多項式(Jacobi Polynomials)\(P_n^{(\alpha,\beta)}(x)\)
・エルミート多項式(Hermite Polynomials):\(H_n(x)\)
とその1階微分、2階微分

のfortran90用のプログラムです。
複素平面で有効です。

このプログラムを使用して生じた問題に対する責任は一切取りません。
ミスがあるかもしれないことを念頭に置いて使ってください。
また、使えない場合がありましたらご連絡いただけると幸いです。

※2017/10/27
整数次数の、実数、複素引数エルミート多項式を追加しました。

特殊関数のモジュール”specialfunction”

プログラム例

呼び出す際は総称名
関数：Laguerre(n,a,x), legendre(l,m,x), jacobi(n,a,b,x)
1階微分：Laguerre1(n,a,x), legendre1(l,m,x), jacobi1(n,a,b,x)
2階微分：Laguerre2(n,a,x), legendre2(l,m,x), jacobi2(n,a,b,x)
を用いてください。

ラゲール陪多項式\(L_n^{(a)}(x)\)の引数の型(n,a,x)の組み合わせは、
(n,a,x) : 整数型、　整数型、　　倍精度型　:　出力は倍精度型
(n,a,x) : 整数型、倍精度型、　　倍精度型　:　出力は倍精度型
(n,a,x) : 整数型、倍精度型、複素倍精度型　:　出力は複素倍精度型

ルジャンドル陪多項式\(P_l^m(x)\)の引数の型(l,m,x)の組み合わせは、
(l,m,x) : 整数型、整数型、　　倍精度型　:　出力は倍精度型
(l,m,x) : 整数型、整数型、複素倍精度型　:　出力は複素倍精度型

ヤコビ多項式\(P^{(a,b)}_n(x)\)の引数の型(n,a,b,x)の組み合わせは、
(n,a,b,x) : 整数型、　　倍精度型、　　倍精度型、複素倍精度型　:　出力は複素倍精度型
(n,a,b,x) : 整数型、複素倍精度型、複素倍精度型、複素倍精度型　:　出力は複素倍精度型

というものが有効になっています。

実際に全部使っているプログラムはこちら。

program main
  use specialfunction
  implicit none
 
  write(6,*)laguerre(4,2,4.78d0)
  write(6,*)laguerre(3,2.5d0,4.78d0)
  write(6,*)laguerre(8,dcmplx(1.6d0,3.01d0),dcmplx(5.64d0,49.1d0))
 
  write(6,*)legendre(8,5,0.28d0)
  write(6,*)legendre(6,2,dcmplx(1.8d0,2.7d0))

  write(6,*)jacobi(3,2.15d0,1.85d0,dcmplx(4d0,8.1d0))
  write(6,*)jacobi(3,dcmplx(1.2d0,0.5d0),dcmplx(5d0,0.3d0),dcmplx(4d0,8.1d0))

  stop
end program main

引数は適当に入れています。
spfunc.f90にmoduleを入れ、
temp.f90に上の12行のコードを書きコンパイル→実行すると、

$ ifort spfunc.f90 temp.f90
$ ./a.out    
    3.2997056066666568    
  -8.4458666666671553E-002
 (  454556598.95966828     ,  376437598.22472727     )
   9378.9803534037801    
 ( -519656.04147187521     ,  135464.19972749997     )
 ( -10936.767937499992     , -2092.2704999999974     )
 ( -18248.193121999979     , -10218.566324666655     )
$

と言う結果を得ます。
wolfram alphaと確かめてみましょう。計算で入れた値はそれぞれ、
\(
\begin{align}
& L_4^2(4.78) \\
& L_3^{2.5}(4.78) \\
& L_8^{1.6+i3.01}(5.64+i49.1) \\
& P_8^{5}(0.28) \\
& P_6^{2}(1.8+i2.7) \\
& P_3^{(2.15,1.85)}(4+i8.1) \\
& P_3^{(1.2+i0.5,5+i0.3)}(4+i8.1) \\
\end{align}
\)
です。

wolframの結果と比較したのがこちら。
spfunc_wolfram
下線部が上記プログラムとwolframとで一致した桁です。
まぁまぁいいんじゃないでしょうか。
\(L_8^{1.6+i3.01}(5.64+i49.1)\)の結果が著しく良くないのが気になります･･･使う際は安定な範囲を確かめてから使ってください。

以下、ここで示している特殊関数の定義と関係式です。

ラゲール陪多項式

  • 定義域

    \(\;\;\;\;
    0 \lt x \lt \infty
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle
    \left[x\frac{d^2}{dx^2}+(\alpha+1-x)\frac{d}{dx}+n\right]L_m^{(\alpha)}(x)=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & L^{(0)}_0(x)=1 \\
    & L^{(0)}_1(x)=-x+1 \\
    & L^{(0)}_2(x)=\frac{1}{2}(x^2-4x+2) \\
    & L^{(0)}_3(x)=\frac{1}{6}(-x^3+9x^2-18x+6) \\
    & L^{(1)}_0(x)=1 \\
    & L^{(1)}_1(x)=-x+2 \\
    & L^{(1)}_2(x)=\frac{1}{2}\left(x^2-6x+6\right) \\
    & L^{(1)}_3(x)=\frac{1}{6}\left(-x^3+12x^2-36x+24\right)
    \end{align}
    \)

  • 漸化式

    \(\;\;\;\;
    \begin{align}
    L_n^{(\alpha)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
    &\hspace{1em}a_{m-1}(x)=1-\frac{n-m+1}{m(\alpha+m)}\cdot x\cdot a_m(x) \\
    &\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
    \end{align}
    \)
    [1]
    1階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d}{dx}L_n^{(\alpha)}(x)=-L_{n-1}^{(\alpha+1)}(x)
    \end{align}
    \)
    2階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d^2}{dx^2}L_n^{(\alpha)}(x)=L_{n-2}^{(\alpha+2)}(x)
    \end{align}
    \)

  • 直交性

    \(\;\;\;\;
    \displaystyle
    \int_{0}^{\infty} e^{-x}x^{\alpha}L^{(\alpha)}_n(x)L^{(\alpha)}_{n’}(x) dx =\frac{\Gamma(\alpha+n+1)}{n!}\delta_{n,n’},\;\;(\alpha>-1)
    \)
    [2]

ルジャンドル陪多項式

  • 定義域

    \(\;\;\;\;
    1 \lt x \lt -1
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle (1-x^2)\frac{d^2}{dx^2}P_l^m(x)-2x\frac{d}{dx}P_l^m(x)+\left[l(l+1)-\frac{m^2}{1-x^2}P_l^m(x)\right]=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & P_0^0(x)=1 \\
    & P_{1}^{-1}(x)=\frac{1}{2}(1-x^2)^{1/2} \\
    & P_{1}^{0}(x)=x \\
    & P_{1}^{1}(x)=-(1-x^2)^{1/2} \\
    & P_{2}^{-2}(x)=\frac{1}{8}(1-x^2) \\
    & P_{2}^{-1}(x)=\frac{1}{2}x(1-x^2)^{1/2} \\
    & P_{2}^{0}(x)=\frac{1}{2}(3x^2-1) \\
    & P_{2}^{1}(x)=-3(1-x^2)^{1/2} \\
    & P_{2}^{2}(x)=3(1-x^2)
    \end{align}
    \)

  • 漸化式

    用いている式
    \(\;\;\;\;
    \begin{align}
    & P_{|m|}^{|m|}(x)=(-1)^{|m|}(2|m|-1)!!(1-x^2)^{\frac{|m|}{2}} \\
    & P_{l}^{-m}(x)=(-1)^m\frac{(l-m)!}{(l+m)!}P_l^m(x) \\
    & P_{|m|}^{|m|+1}(x)=(2|m|+1)xP_{|m|}^{|m|}(x) \\
    & P_{|m|}^{|m|+q}(x)=\left(\frac{2|m|-1}{q}+2\right)xP_{|m|}^{|m|+q-1}(x)-\left(\frac{2|m|-1}{q}+1\right)P_{|m|}^{|m|+q-2}(x)
    \end{align}
    \)[3]
    1階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1&\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=\frac{lx}{(x^2-1)}P_l^m(x)-\frac{l+m}{(x^2-1)} P_{l-1}^m(x) \\
    m\ne\pm 1\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=c_{2}P^{m+2}_{l-1}(x)+c_0 P^{m}_{l-1}(x)+c_{-2} P^{m-2}_{l-1}(x) \\
    &\hspace{1em} c_{2}=\frac{1}{4(m+1)}\\
    &\hspace{1em} c_{0}=\frac{l+m}{2}\left(1+\frac{l}{1-m^2}\right)\\
    &\hspace{1em} c_{-2}=-\frac{(l+m)(l+m-1)(l+m-2)(l-m+1)}{4(m-1)}
    \end{align}
    \)
         ※\(m=\pm 1\)の時、\(x=\pm 1\)で発散します。
    2階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=\frac{1}{1-x^2}\left[2x\frac{dP_l^m(x)}{dx}+\frac{(l+1)(l+m)}{2l+1}\frac{dP_{l-1}^m(x)}{dx}-\frac{l(l-m+1)}{2l+1}\frac{dP_{l+1}^m(x)}{dx}\right] \\
    m\ne\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=c_{2}\frac{dP^{m+2}_{l-1}(x)}{dx}+c_0 \frac{dP^{m}_{l-1}(x)}{dx}+c_{-2} \frac{dP^{m-2}_{l-1}(x)}{dx}
    \end{align}
    \)
         ※\(m=\pm 1,\pm 3\)の時、\(x=\pm 1\)で発散します。

  • 直交性

    \(l\)に対する直交性
    \(\;\;\;\;
    \displaystyle
    \int_{-1}^{1} P_{m}^{l}(x)P_{m}^{l’}(x) dx =\frac{2}{2l+1}\cdot \frac{(l+m)!}{(l-m)!}\delta_{l,l’}
    \)
    [2]

ヤコビ多項式

  • 定義域

    \(\;\;\;\;
    -1 \lt x \lt 1
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle
    \left[(1-x^2)\frac{d^2}{dx^2}+(\beta-\alpha-(\alpha+\beta+2)x)\frac{d}{dx}+n(n+\alpha+\beta+1)\right]P_n^{(\alpha,\beta)}(x)=0
    \)

  • 漸化式

    \(\;\;\;\;
    \begin{align}
    P_n^{(\alpha,\beta)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
    &\hspace{1em}a_{m-1}(x)=1-\frac{(n-m+1)(\alpha+\beta+n+m)}{2m(\alpha+m)}(1-x)\cdot a_m(x)\\
    &\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
    \end{align}
    \)
    [1]
    1階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d}{dx}P_n^{(\alpha,\beta)}(x)=\frac{\Gamma(\alpha+\beta+n+2)}{2\Gamma(\alpha+\beta+n+1)}P_{n-1}^{(\alpha+1,\beta+1)}(x)
    \end{align}
    \)
    2階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d^2}{dx^2}P_n^{(\alpha,\beta)}(x)=\frac{\Gamma(\alpha+\beta+n+3)}{2^2\Gamma(\alpha+\beta+n+1)}P_{n-2}^{(\alpha+2,\beta+2)}(x)
    \end{align}
    \)

  • 直交性

    \(\;\;\;\;
    \begin{align}
    \displaystyle
    &\int_{-1}^{1} (1-x)^{\alpha}(1+x)^{\beta}P_{n}^{(\alpha,\beta)}(x)P_{n’}^{(\alpha,\beta)}(x)dx \\
    &\hspace{4em}=\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}
    \frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{n!\Gamma(n+\alpha+\beta+1)}\delta_{nm}\\
    &(\alpha>-1,\ \beta>-1)
    \end{align}
    \)    [2]

参考文献

[1] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.18, Page 789)
[2] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.2, Page 774,775)
[3] Associated Legendre polynomials -wikipedia

数学

平方根を手計算で簡単に求める方法

2件のコメント

平方根の値を簡単に得る方法です。

もしかしたら、某福〇漫画に出てくるように、\(\sqrt{2}\)を計算しないと命がヤバいみたいな状況に陥るかもしれません。そんな状況を乗り切りましょう。

開平法による筆算とかニュートン法を手でやる方法とかありますが、
覚えるのが大変ですし、計算過程が大変でミスが起こります。

使うべきは
連分数
です。
連分数を使うと平方根の計算が1回の割り算に変わります(確かめのためには2回の割り算が必要)。

早速本題に入りましょう。
\(\;\;\;\;\sqrt{2}\)
を求めてみます。

まず式変形をします。
\(
\;\;\;\;\sqrt{2}=1+(\sqrt{2}-1)
\)
として、整数部分\(1\)と小数部分\((\sqrt{2}-1)\)に分けます。
そして、小数部分に対し、
\(
\;\;\;\;\displaystyle \frac{\sqrt{2}+1}{\sqrt{2}+1}
\)
を掛けますと、
\(
\;\;\;\;\displaystyle \sqrt{2}=1+\frac{1}{1+\sqrt{2}}
\)
になります。ここで、整数部分を左辺に移行しまして、
\(
\;\;\;\;\displaystyle \sqrt{2}-1=\frac{1}{1+\sqrt{2}}
\)
とします。ここで、
\(
\;\;\;\;x=\sqrt{2}-1
\)
と置いて、\(\sqrt{2}=x+1\)を代入すると
\(\;\;\;\;
\begin{align}
x&=\frac{1}{2+x} \\
x&=(2+x)^{-1} \;\;\cdots (1)
\end{align}
\)
なる関係式が得られます。
\(\sqrt{2}=x+1\)ですから、\(x\)が求まればいいわけです。

ここまで来れば後は簡単です。漸化式と見立てて(1)を逐次的に解けばいいのです。
まず、(1)式に\(x(=\sqrt{2}-1)\)の近似値として0という値を入れます。
この値は\(x\)の近似値であれば何でも構いませんが、後々のために整数か分数にしたほうがいいです。

(1)の右辺に代入します。{　}の中身は、そこで打ち切った時に割り算して求めた場合の結果で,アンダーラインは\(x=\sqrt{2}-1=0.414213562372\cdots\)との一致具合を表します。

[1回目]
\(
\displaystyle (2+0)^{-1}=\frac{1}{2}\;\;\;\{\cdots x=0.5\}
\)
[2回目]
\(
\displaystyle (2+\frac{1}{2})^{-1}=\frac{2}{5}\;\;\;\{: x=0.\underline{4}\}
\)
[3回目]
\(
\displaystyle (2+\frac{2}{5})^{-1}=\frac{5}{12}\;\;\;\{: x=0.\underline{41}66\cdots\}
\)
[4回目]
\(
\displaystyle (2+\frac{5}{12})^{-1}=\frac{12}{29}\;\;\;\{: x=0.\underline{41}379\cdots\}
\)
[5回目]
\(
\displaystyle (2+\frac{12}{29})^{-1}=\frac{29}{70}\;\;\;\{: x=0.\underline{4142}857\cdots\}
\)
[6回目]
\(
\displaystyle (2+\frac{29}{70})^{-1}=\frac{70}{169}\;\;\;\{: x=0.\underline{4142}011\cdots\}
\)
[7回目]
\(
\displaystyle (2+\frac{70}{169})^{-1}=\frac{169}{408}\;\;\;\{: x=0.\underline{41421}568\cdots\}
\)
[8回目]
\(
\displaystyle (2+\frac{169}{408})^{-1}=\frac{408}{985}\;\;\;\{: x=0.\underline{414213}131\cdots\}
\)
[9回目]
\(
\displaystyle (2+\frac{408}{985})^{-1}=\frac{985}{2378}\;\;\;\{: x=0.\underline{414213}624\cdots\}
\)
[10回目]
\(
\displaystyle (2+\frac{985}{2378})^{-1}=\frac{2378}{5741}\;\;\;\{: x=0.\underline{4142135}551\cdots\}
\)
[11回目]

\(
\displaystyle (2+\frac{2378}{5741})^{-1}=\frac{5741}{13860}\;\;\;\{: x=0.\underline{41421356}421\cdots\}
\)
…

ただ2倍して足して、分子分母をひっくり返せば計算できるので、何も難しいことはありません。
ここまでの計算は5分もかからないでしょう。

しかも、途中で計算ミスをしても構わないのです。
なぜなら、近似値からさらに良い値にどんどん更新していくからです。計算回数を多くしないといけないだけです。

収束性を確かめたいのなら、もう一回増やしたものを計算する必要があります。その差を見ると良いでしょう。

一回だけ割り算をしなければなりませんね。これはどうしてもなくせません。
ただし、\(\sqrt{2}\)が1.41421356と分かっていれば、筆算の計算は途中までは掛け算と引き算で行けます。そう考えるとまぁ多少は簡単なんですかね・・・？

参考文献

連分数の応用 -講義ノート、筑波大

gnuplot, 数学

正多角形とスピログラフの数式

6件のコメント

正多角形とスピログラフの数式です。

  1. 正多角形(2次元)
  2. 正多角形(3次元)
  3. 星型多角形(2次元)
  4. 星型多角形(3次元)
  5. トロコイド(2次元)
  6. トロコイド(3次元)
  7. スピログラフ

2016/10/24
幾何学模様を作りたい方は幾何学模様(gnuplot)へどうぞ。

[adsense1]

正多角形の式(2次元)

正n角形の数式は[1]を参考に修正しました。
正n角形の数式は、
\(
\begin{equation}
\left\{
\begin{aligned}
x&\displaystyle =\cos t \cdot \frac{\cos \frac{a_n}{2}}{\cos\left[a_n\left\{\frac{t}{a_n}-{\rm floor}\left(\frac{t}{a_n}\right)\right\}-\frac{a_n}{2}\right]}\\
y&\displaystyle =\sin t \cdot \frac{\cos \frac{a_n}{2}}{\cos\left[a_n\left\{\frac{t}{a_n}-{\rm floor}\left(\frac{t}{a_n}\right)\right\}-\frac{a_n}{2}\right]}
\end{aligned}
\right.
\end{equation}
\)
です。ここで、\(a_n=\frac{2\pi}{n}\)と置きました。また、\({\rm floor}(x)\)は床関数(床関数と天井関数)を表します。
媒介変数\(t\)の範囲は(\(0\le t \lt 2\pi\))です。
グラフを書かせてみますとこうなります。

polygon3_10_c
n=3,4,5,6
n=7,8,9,10の順で、単位円を一緒に記述しています。
線の色の濃さは媒介変数の値に応じて変化させています。

gnuplotのスクリプトはこちら

多面体の式(3次元)

半径1の多面体は、平面のx,yの媒介変数による正n角形の関数をそれぞれpc(t),ps(t)と置くと、
\(
\begin{equation}
\left\{
\begin{aligned}
x& = {\rm ps}(u)\cdot {\rm pc}(v) \\
y& = {\rm ps}(u)\cdot {\rm ps}(v) \\
z& = {\rm pc}(u)
\end{aligned}
\right.
\end{equation}
\)
と書けます。
媒介変数\(u,v\)の範囲は(\(0\le u \lt \pi, \;\; 0\le v \lt 2\pi\))です。

プロットすると以下のようになります。
polygon3d_3_10_c
n=3,4,5,6
n=7,8,9,10の順です。

gnuplotスクリプトはこちら

星型多角形(2次元)

\(
\begin{equation}
\left\{
\begin{aligned}
x&\displaystyle =\cos t \cdot \frac{\cos a_n}{\cos\left[a_n -\frac{1}{n}\cos^{-1}\{\cos(nt)\}\right]}\\
y&\displaystyle =\sin t \cdot \frac{\cos a_n}{\cos\left[a_n -\frac{1}{n}\cos^{-1}\{\cos(nt)\}\right]}
\end{aligned}
\right.
\end{equation}
\)
です(\(n\ge 5\))。ここで、\(a_n=\frac{2\pi}{n}\)と置きました。
グラフを書かせてみるとこのようになります。
starpolygonp_c

gnuplotスクリプトはこちら。
※\(\cos^{-1}(\cos(x))\)は不安定です。
場合により書けなくなることがあります。ここでは上記関数の代わりに、正弦波ではない周期関数の数式
で記述した関数\(s(a,x)\)を用いています。

星型多角形(3次元)

3次元の場合は多面体の時と同じように、平面のx,yの媒介変数による星型多角形の関数をそれぞれstpc(t),stps(t)と置くと、
\(
\begin{equation}
\left\{
\begin{aligned}
x& = {\rm stps}(u)\cdot {\rm stpc}(v) \\
y& = {\rm stps}(u)\cdot {\rm stps}(v) \\
z& = {\rm stpc}(u)
\end{aligned}
\right.
\end{equation}
\)
と書けます。
媒介変数\(u,v\)の範囲は(\(0\le u \lt \pi, \;\; 0\le v \lt 2\pi\))です。
グラフはこのようになります。
starpolygon3dp_c

gnuplotスクリプトはこちら。

トロコイド(2次元)

トロコイドと呼ばれる曲線があります[2]。
これをgnuplotで書いてみましょう。数式は[2]を参考にして以下の通りです。
\(
\begin{equation}
\left\{
\begin{aligned}
x& = \left[(n-1)\cos t +h\cos((n-1)t)\right]/n \\
y& = \left[(n-1)\sin t +h\sin((n-1)t)\right]/n
\end{aligned}
\right.
\end{equation}
\)
\(h=1\)の時、半径が1の円に内接するように決めています。
また、\(n\)が整数の時、内接する円が外側の円をちょうど一周回った時に元の位置に戻る条件です。
別に整数でなくても構いません。その時はもっと複雑なものになることでしょう。

書いてみますと曲線はこうなります。
自由度\(h\)を1にしてグラフを書いてみましょう。
trocoid2d_3_10_c

自由度hを少し変えてみましょう。
trocoid2dh_3_10_c

トロコイド(3次元)

3次元の場合は多面体の時と同じように、平面のx,yの媒介変数によるトロコイドの関数をそれぞれhcyc(t),hcys(t)と置くと、
\(
\begin{equation}
\left\{
\begin{aligned}
x& = {\rm hcys}(u)\cdot {\rm hcyc}(v) \\
y& = {\rm hcys}(u)\cdot {\rm hcys}(v) \\
z& = {\rm hcyc}(u)
\end{aligned}
\right.
\end{equation}
\)
と書けます。
媒介変数\(u,v\)の範囲は(\(0\le u \lt \pi, \;\; 0\le v \lt 2\pi\))です。

trocoid3d_3_10_c
n=3,4,5,6
n=7,8,9,10の順です。

gnuplotスクリプトはこちら

[adsense2]

スピログラフ(2次元)

スピログラフは遊んだ方も多いと思います。これです。
1024px-Spirograph
©Alexei Kouprianov/2007/CC-BY 2.5
スピログラフはトロコイドの特別な時です。スピログラフをもっと一般化したのがトロコイドです[2]。

図示すれば、スピログラフと言うのは
speq_compressed
という二つの円からなる軌道です。黄色の点線の長さが等しい、という束縛条件が課せられます。
数式は上のトロコイドと同じで
\(
\begin{equation}
\left\{
\begin{aligned}
x& = \left[(n-1)\cos t +h\cos((n-1)t)\right]/n \\
y& = \left[(n-1)\sin t +h\sin((n-1)t)\right]/n
\end{aligned}
\right.
\end{equation}
\)
と書けます。

もう少し考えを一般化して、周りを取り囲んでいる形状が円ではない時を考えましょう。その時の数式を考えます。
考え方はトロコイドと一緒ですが、もう少し拡張しましょう。

今、2つの媒介変数による軌道を持っていたとします。１つの媒介変数の軌道上を動かせばスピログラフの考えと同じになります。図で示すと以下のようになります。
spirograph_any_c
スピログラフは両方円であり、滑りが無い(cを決定する条件)場合に対応するわけです。ただ、円である必要性はないので、多角形であろうが、何だろうが数式上では書けるわけです。
興味深い事に、大きな軌道に沿って小さな軌道を動かした、と見ようが小さな軌道に沿って大きな軌道を動かしたとみても数式上では一緒なのです。面白いですね。

いくつか例を挙げます。
五角形上を動かす場合はこうです。
parameter1

数式として,
\(
\begin{equation}
\left\{
\begin{aligned}
x& = f_{1x}(t)+r\cdot f_{2x}(s\cdot t) \\
y& = f_{1y}(t)+r\cdot f_{2y}(s\cdot t)
\end{aligned}
\right.
\end{equation}
\)
と書き表すことにします。これを元にgnuplotのスクリプトを書けば以下のようになります。

色々に変化させた時の色々なときのグラフです。
spr2_any_compressed

rを綺麗な数字にしなければもっと複雑な図が得られます。
また、この考え方は3つの媒介変数の軌道、4つの媒介変数の軌道…として容易に拡張できます。
もっと複雑な模様が欲しければ、もう少し増やしてみましょう。

例えば3つの場合はこうなります。
parameter3_spi
このスクリプトを置いておきます。

この考えを拡張していくと”フーリエ級数の図示化”と言うのができるわけです。

実際にやっている方がいるので紹介しておきますね。

参考文献

[1]Is there an equation to describe regular polygons? -MATHEMATICS
[2]Hypotrochoid -Wolfram Mathworld

[3]Is there an equation to describe regular polygons? -Mathematics
2021/11/29追記）滑らかに星と円を遷移する場合、このページが丁寧です。

gnuplot, 数学

アニメーションに便利な関数

コメントする

アニメーションに便利な関数、手法です。

点の表示

plot sprintf("< echo ''%f %f''", x0,y0)

とすると点\((x_0,y_0)\)が描写されます。

窓関数

窓関数です。関数\(W(x)\)は、
\(
\displaystyle W(x)=\frac{r^{2^n}}{(x-x_c)^{2^n}+r^{2^n}}
\)
書かれます。この関数は\(r\)が窓の横幅を示し、\(n\)が窓の端の様子を決めるパラメータです。窓の広がりはおよそ\(2r\)です。

fwi(n,r,i,ic)=(real(r)**(2**n))/((real(i-ic))**(2**n)+real(r)**(2**n))

加速、減速

extanh
\(
\begin{align}
y(x)&=C_0\tanh(a(x-x_c))+C_1 \\
& C_0 = \frac{f_0-f_1}{t_0-t_1} \\
& C_1 = \frac{-t_1f_0+t_0f_1}{t_0-t_1} \\
& t_0 = \tanh(a(x_0-x_c)) \\
& t_1 = \tanh(a(x_1-x_c)) \\
& x_c = \frac{x_0+x_1}{2}
\end{align}
\)
の関数は点\((x_0,y_0)\)から点\((x_1,y_1)\)に\(\tanh(ax)\)で与えられる関数で書きます。
関数の中心は\((x_c,y_c)=((x_0+x_1)/2,(y_0+y_1)/2)\)になります。
画像の2つの赤い点は\((x_0,y_0)\)と\((x_1,y_1)\)を表しています。
\(a\)は変化の度合いを調節するパラメータです。

この関数は非常に使えます。ラミエルを作っている最中にこの関数を考えましたが、想像以上に役立ちます。

fth(a,x,x0,f0,x1,f1)=(f0-f1)*tanh(a*(real(x)-(x0+x1)*0.5e0))/(tanh(a*(real(x0)-(x0+x1)*0.5e0))-tanh(a*(real(x1)-(x0+x1)*0.5e0)))+(-f0*tanh(a*(real(x1)-(x0+x1)*0.5e0))+f1*tanh(a*(real(x0)-(x0+x1)*0.5e0)))/(tanh(a*(real(x0)-(x0+x1)*0.5e0))-tanh(a*(real(x1)-(x0+x1)*0.5e0)))
 # x0,f0 --> x1,f1

この関数を用いた応用です。
extime_tanh

これはある時間内で高さを反転させるものです。コードはこちら。

set xr[-6:6]
set yr[0.8:2.2]
set size ratio 0.5 1
   
a=1e0
at=0.4e0
x0=-5e0
f0=1e0
x1=5e0
f1=2e0

#set term gif enhanced animate optimize delay 8 size 800,500
#set output 'extime_tanh.gif'

i0=0
i1=20
do for[i=i0:i1]{
    ft0=fth(at,i,i0,f0,i1,f1)
    ft1=fth(at,i,i0,f1,i1,f0)
    plot fth(a,x,x0,ft0,x1,ft1) w l lt 2 lc 3 lw 5 ti sprintf("a=%3.2f, at=%3.2f, (x_0,f_0)=(%3.2f,%3.2f), (x_1,f_1)=(%3.2f,%3.2f)",a,at,x0,f0,x1,f1), sprintf("< echo ''%f %f''", x0,ft0) pt 7 ps 2 lc 7, sprintf("< echo ''%f %f''", x1,ft1) pt 7 ps 2 lc 7

  }

#set out
#set terminal wxt enhanced