「プログラミングと数値計算」カテゴリーアーカイブ

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

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

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

テキストエディタ, プログラミングと数値計算

emacs上で再読み込み

コメントする

emacs上でファイルを再読み込みしたいとします。
OSはLinux mint17.1です。

例えば、
emacsであるファイルを開いたまま、sedコマンドでファイルをいじると、そのいじったファイルは現在開いているemacsには反映されません。

これを変更してF5キーを押すだけでemacs上で再読み込みができるようにします。

ホームディレクトリに
.emacs
というファイルを作成し、その中に↓

(defun revert-buffer-no-confirm (&optional force-reverting)
  "Interactive call to revert-buffer. Ignoring the auto-save
 file and not requesting for confirmation. When the current buffer
 is modified, the command refuses to revert it, unless you specify
 the optional argument: force-reverting to true."
  (interactive "P")
  ;;(message "force-reverting value is %s" force-reverting)
  (if (or force-reverting (not (buffer-modified-p)))
      (revert-buffer :ignore-auto :noconfirm)
    (error "The buffer has been modified")))

  ;; reload buffer
  (global-set-key (kbd "<f5>") 'revert-buffer-no-confirm)

と入力します。これを記述し、もう一度emacsを立ち上げ,sedコマンド等で外部から編集し、F5キーを押せば更新されます。

参考文献

Emacsで現在開いてるバッファを再読込する

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

sedコマンドによるファイルの書き換え

コメントする

シェルスクリプトによってファイルを書き換えます。

  1. sedコマンドの概要
  2. 具体例一覧
    1. 特定の数値を変えたい場合
    2. 特殊文字(/等)を含む文字列の変更
    3. 置換する数値を変数での指定
    4. ある数値以降を全て変更したい
    5. ファイルの空白行の削除
    6. ファイルの行頭に数値を入れる
    7. ファイルの最後に改行を入れる

      8. 繰り返し変更したい場合

    8. do文で何回も変更したい
    9. 全く違う変数を繰り返し代入したい
    10. 整数の変数を規則的に何回も変更したい
    11. 小数点を含む変数を規則的に何回も変更したい
    12. あるディレクトリ以下にある全ファイルを変数にしたい

sedコマンドの概要

対象とするファイル名”input.ini”の中身に

&input
 a=7,
 b=13,
 c='./dat/',
&end

という文が書いてあるとします。
これのファイルの中身を端末からのコマンドで書き換えるためにはsedコマンドを利用します。

sedコマンドとは、置換用のコマンドです。
例えば、コマンドライン上で

sed -e "s/a=7/a=9/" ./input.ini

と入力すると”input.ini”の中身の”a=7″に一致する場所を”a=9″に置換して出力します。
ここで重要なのが、”a=7″に一致する場所だけなので、後ろにあるコンマは影響されません。
実際に実行してみると

$ sed -e "s/a=7/a=9/" ./input.ini
&input
 a=9,
 b=13,
 c='./dat/',
&end

となってコマンドライン上に出力されます。
この時、書き換わった内容はコマンドラインのみに出力されるため、元のファイル自体に書き換えはされません。

元のファイルの書き換えを行いたければ、出力先を適当なファイル、例えば”tmp”というファイルにし、それを”input.ini”にリネームすればいいです。
すなわち、

sed -e "s/a=7/a=9/" ./input.ini > tmp
mv tmp input.ini

とすればいいのです。実際に実行して中身を見ますと

$ sed -e "s/a=7/a=9/" ./input.ini > tmp
$ mv tmp input.ini
$ cat input.ini
&input
 a=9,
 b=13,
 c='./dat/',
&end

となります。

sedコマンドの概要はここで終了です。
以降は具体例とシェルスクリプトのコード(ファイル名”q.sh”),実行例を載せていきます。

特定の数値を変えたい場合

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s/a=7/a=9/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=9,
 b=13,
 c='./dat/',
&end

特殊文字(/等)を含む文字列の変更

sedコマンドを使うために/を使います。ならば/を置換したい場合どうすればいいでしょう。
この時はsedコマンドで使っている区切りを表す/を別の特殊文字(例えば%)に置き換えればよいです。

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s%c='./dat/'%c='./data/'%" input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=7,
 b=13,
 c='./data/',
&end

置換する数値の変数での指定

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
pr=5
sed -e "s/a=7/a=${pr}/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=5,
 b=13,
 c='./dat/',
&ends

ある数値以降を全て変更したい

コンマも含めて消したいとします。
対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s/c=.*/c=12345/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=7,
 b=13,
 c=12345
&end

ファイルの空白行の削除

ファイルに含まれる空白行を全て削除します。
対象とするファイル”input.ini”

&input
 a=7,

 b=13,

 c='./dat/',


&end

“q.sh”の中身

#!/bin/sh
sed -e '/^ *$/d' input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=7,
 b=13,
 c='./dat/',
&end

ファイルの行頭に数値を入れる

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 1000 ]
do
    sed -e "${i}s/^/${i} /" input.ini > temp
    mv temp input.ini

    i=`expr $i + 1`
done

実行例

$ sh q.sh
$ cat input.ini
1 &input
2  a=7,
3  b=13,
4  c='./dat/',
5 &end
6

“q.sh”中、1000をファイルに応じて変更してください。

※また、行番号を入れるだけの場合、catコマンドを利用して、

cat -b input.ini > tmp
mv tmp input.ini

とするのが簡単だと思います。
オプション”-b”は空行を含めず行番号を付けるという意味です。空行を含めたい場合は”-n”でできます。

ファイルの最後に改行を入れる

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed '$a\\n' input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
 a=7,
 b=13,
 c='./dat/',
&end
   

$

※どうしても2行分改行されてしまうようです。
sedコマンドで1行分の改行は、どうすればよいか分かりませんでした。

苦肉の策ですが、適当なファイル(例えばnn.txt)に改行だけを入力し、catコマンドでファイルを結合すると一応できます。
—nn.txt—

————
を作っておいて、
#!/bin/sh
cat input.ini nn.txt > tmp
mv tmp input.ini
とすれば”input.ini”の後ろに”nn.txt”が結合されるため、1行分の改行が可能です。

do文で何回も変更したい

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
for pr in 0 1 2
do
    sed -e "s/a=.*/a=${pr},/" input.ini > tmp
    mv tmp input.ini
   
    cat input.ini
done

実行例

$ sh q.sh
&input
 a=0,
 b=13,
 c='./dat/',
&end
&input
 a=1,
 b=13,
 c='./dat/',
&end
&input
 a=2,
 b=13,
 c='./dat/',
&end

全く違う変数を繰り返し代入したい

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
for j in 3d0 939d0
do
    echo "count -> ${i}"
    sed -e "s/a=.*/a=${j},/" ./input.ini > tmp
    mv tmp input.ini
   
    grep "a=" input.ini
   
    i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
 a=3d0,
count -> 2
 a=939d0,
$

シェルスクリプトの中の”grep”は”a=”に一致する行を書き出させるものです。

整数の変数を規則的に何回も変更したい

対象とするファイル”input.ini”

&input
 a=7,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 3 ]
do
    echo "count -> ${i}"
    sed -e "s/a=.*/a=${i},/" ./input.ini > tmp
    mv tmp input.ini
   
    grep "a=" input.ini
   
    i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
 a=1,
count -> 2
 a=2,
count -> 3
 a=3,
$

小数点を含む変数を規則的に何回も変更したい

小数点を含む演算は”expr”コマンドではできません。
シェルスクリプトで数値計算を行う”bc”コマンドを用いましょう。
ここでは、変数\(v\)に、
\(v=i\times0.7,\;\;(i=1,2,3)\)
を代入する演算を考えましょう。

対象とするファイル”input.ini”

&input
 a=3,
 b=13,
 c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 3 ]
do
    echo "count -> ${i}"

    v=$(echo "scale=3; $i*0.70 " | bc | sed -e 's/^\./0./g')
    echo "$v"
   
    sed -e "s/a=.*/a=${v},/" ./input.ini > tmp
    mv tmp input.ini
   
    grep "a=" input.ini
   
    i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
0.70
 a=0.70,
count -> 2
1.40
 a=1.40,
count -> 3
2.10
 a=2.10,
$

あるディレクトリ以下にある全ファイル名を変数にしたい

状況

.
├── a.sh
├── dat
├── data
│   ├── asdfgh.d
│   ├── qwerty.d
│   └── zxcvbn.d
└── input

inputの中身

&input
  filename="./data/xxxxx.d",
  key = 0,
&end

a.shの中身

#!/bin/sh
i=1
for j in ./data/*
do
    sed -e "s%filename=.*%filename="${j}",%" input > tmp
    mv tmp ./dat/${j##*/}_ex
    echo ${j##*/}
    grep "filename=" ./dat/${j##*/}_ex
done

———————-

a.sh実行

$sh a.sh
asdfgh.d
  filename="./data/asdfgh.d",
qwerty.d
  filename="./data/qwerty.d",
zxcvbn.d
  filename="./data/zxcvbn.d",

a.sh実行後の構成

.
├── a.sh
├── dat
│   ├── asdfgh.d_ex
│   ├── qwerty.d_ex
│   └── zxcvbn.d_ex
├── data
│   ├── asdfgh.d
│   ├── qwerty.d
│   └── zxcvbn.d
└── input

例えばdat/asdfgh.d_exの中身

&input
  filename="./data/asdfgh.d",
  key = 0,
&end

参考先

フィルタを使用した文字列操作 1 -UNIX & Linux コマンド・シェルスクリプト リファレンス
sed の使い方 -サーバエンジニアの知恵袋
制御構文 -シェルスクリプト入門
How to pass results of bc to a variable – Ask Ubuntu
bcによる少数の演算 -SWEng TIPs

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

openMPによる並列計算

コメントする

まとめ

fortran90では、

program main
  implicit none
  integer::i,omp_get_thread_num

  ! number of threads for parallel computation
  call omp_set_num_threads(2)

  !$omp parallel do
  do i=1,100
     write(6,*) omp_get_thread_num()
  enddo
  !$omp end parallel do  

  stop
end program main

とすれば2個のスレッドで並列計算され、コンパイル時に
intel® fortranでは、

ifort -openmp main.f90

gfortranでは、

gfortran -fopenmp main.f90

とすれば良いです。

fortran90でのお話です。

openMPは複数のCPUを使って並列計算をするインターフェースのことです
C/C++,fortranで対応しています。
何も指定しない場合、計算は基本的に1つのコアのみを使います。
Linuxを使っている場合、恐らくデフォルトで使えるようになっているはずです。

並列計算を行うまでの大まかな手順は2つ。

  1. プログラム内部に文章を書く
  2. コンパイル時にオプションを追加
    • intel®fortranコンパイラの場合
    • gfortranコンパイラの場合

です。

※今使っているパソコンの並列計算できるCPUの数の確認は、

cat /proc/cpuinfo | grep processor

です。proc/cpuinfoの場所は使っているOSによって違う可能性があります。
”linux mint cpu 確認”
等とgoogle等で検索するのが手っ取り早いでしょう。

プログラム内部に文章を書く

プログラム内でopenMPを使用するdoループを挟むように、

!$omp parallel do
do i=1,100
    call test(i)   
enddo
!$omp end parallel do

と書きます。基本はこれだけです。
並列計算に使えるスレッドが2個である場合、並列計算というのは

スレッド0

do i=1,50
    call test(i)   
enddo

スレッド1

do i=51,100
    call test(i)   
enddo

と割り振る操作なのです。

上記のようにdoループを分割するため、注意しなければならない点として、
doループに用いる変数に対して別々に計算しても計算結果が変わらない事
が挙げられます。
上記例では、i=0から始めて、i=1,2,3,･･･と順次計算する必要は無く、i=5,10,1,3,･･･みたいな乱雑な順番でも計算結果が変わらない状況で使わなければなりません。
特に、i=2を計算したい時にi=1の計算結果を使う場合は並列化することは出来ません。

実際には同期したりとかで、出来なくはないですが面倒くさく、僕はやったことが無いので書きません。

もしも並列計算したいスレッド数を制限したい場合、

!$omp parallel do

を

!$omp parallel do num_threads(2)

のように指定するか(この例ではその部分の並列計算を2つのスレッドで行う)、
並列計算が始まる前のプログラムの初めの方に

call omp_set_num_threads(2)

と書きましょう。これを書くと全ての並列計算箇所が2個のスレッドで行われることになります。

また、関数omp_get_thread_num()はスレッドにつけられた番号を返す事が出来ます。引数は要りません。
なので、これを利用すると何個のスレッドが使われているのか直に確認できます。

補足ですが、MKL(マス・カーネル・ライブラリ)のルーチン内で使われる並列計算のスレッド数は

call mkl_set_num_threads(2)

と書くことで制御できます。

コンパイル時のオプション

intel® fortranでは、

ifort -openmp main.f90

gfortranでは、

gfortran -fopenmp main.f90

としてコンパイル、実行しましょう。

実例

例として

program main
  implicit none
  integer::i,omp_get_thread_num

  ! number of threads for parallel computation
  call omp_set_num_threads(2)

  !$omp parallel do
  do i=1,4
     write(6,*) omp_get_thread_num()
  enddo
  !$omp end parallel do

  stop
end program main

を動かすと同出力されるか見てみましょう。
gfortranコンパイラで、使えるスレッドの数は2つです。
上記のコードを

gfortran -fopenmp main.f90

としてコンパイルし、実行すると

./a.out
           1
           0
           0
           1

という結果が得られます。
ここに出ている数字はスレッドの番号で、0,1と番号付けされたスレッドが確かに出力されている事が分かります。1,0,0,1の順番はどうでも良いです。実行するごとに変わるので、これは同時並行処理されているためでしょう。

call omp_set_num_threads(1)

として計算を行うと、

./a.out
           0
           0
           0
           0

と表示されることから、確かに使うスレッド数が制御されていることが分かります。

doループの並列処理を処理が終わった順に行う場合

doループ内で、あるときだけ計算が重くなる場合があります。
この場合、均等にdoループの変数を分けてしまうとCPU1つあたりの処理が大きく異なってしまいます。そこで、doループの変数を終わった順に処理していく指定方法があります。それがschedule構文です。以下のようにすると、i=1から計算が終わった順にiを増やしていきます。

program main
  implicit none
  integer::i,omp_get_thread_num,N

  ! number of threads for parallel computation
  call omp_set_num_threads(2)

  N=100
  !$omp parallel do schedule(dynamic,1)
  do i=1,N
     write(6,*) omp_get_thread_num()
  enddo
  !$omp end parallel do
 
  stop
end program main

とします。ただし、1回あたり計算時間がほとんどかからない場合、オーバーヘッド時間が無視できなくなるため、使うときは吟味してから使用するべきです。

また、並列化しても共通の変数を使うことができます。以下のようにすると、計算の進み具合がわかります。

program main
  implicit none
  integer::i,omp_get_thread_num,N,prog

  ! number of threads for parallel computation
  call omp_set_num_threads(2)

  N=100
  prog=0
  !$omp parallel do shared(prog) schedule(dynamic,1)
  do i=1,N
     write(6,*) omp_get_thread_num()
     prog=prog+1
     if(mod(prog,50).eq.0)write(6,'(A,i6,A,i6)')" ",prog," / ",N
  enddo
  !$omp end parallel do
 
  stop
end program main

とします。

経験的に、openmpを使用する際はただ一つのサブルーチンを用いて、並列個所を簡単にするのが良いです。

  !$omp parallel do
  do i=1,N
     a(i)=0
     do j=i,N
         a(i)=a(i)+1
     enddo
     a(i)=a(i)*i     
  enddo
  !$omp end parallel do

と、!$~!$につらつらと書くのではなく、

  !$omp parallel do
  do i=1,N
    call test(i,N,a(i))
  enddo
  !$omp end parallel do

として、サブルーチン１つにするだけの方がバグが出にくいです。

openMPの機能はこれだけではありません。
調べることをお勧めします。ここではこれ以上書きません。

参考文献

インテル® コンパイラーOpenMP*入門

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

MATLABを用いた3次元等値面の描き方

コメントする

さて, 数値計算で得られた計算データを3次元で美しく描きたくなると思います.
なると思います.
MATLABで等値面を描くと, とても綺麗になるのでおすすめです.
ということで, 描き方について紹介したいと思います.
まぁ, 手始めに楕円体ですかね.

program main
  implicit none
 
  integer(8)::na,nb,nc,intangle
  real(8)::x,y,z,V,theta
  integer(8),parameter::no=50
  real(8)::a,b,c
  real(8)::ss(1:3)

  Open(50,FILE='data.dat',STATUS='unknown')
     
     a=0.6d0
     b=1.2d0
     c=0.6d0
     
     do nc=-no,no
        write(*,*)nc+no,"/",no*2
        do na=-no,no
           do nb=-no,no
             
              ss(1)=dble(na)/dble(no)
              ss(2)=dble(nb)/dble(no)
              ss(3)=dble(nc)/dble(no)
             
              V=(b**2d0 *c**2d0 *SS(1)**2d0)+(c**2d0 *a**2d0 *SS(2)**2d0)+(a**2d0 *b**2d0 *SS(3)**2d0) !楕円体の式
             
              write(50,'(3f10.5,2f11.7)')ss(1),ss(2),ss(3),V
           end do
        end do
     end do
     
     close(50)
     
  end do
end program main

今回は, 101x101x101で計算しますが, これらの範囲は自由に変更できます.
実際に使う際は, vの式を任意の式に変更するだけです.

x, y, zの回す順番だけは, このプログラムコードに合わせて下さい.

さて, data.datができました. 次に, MATLABで読み込めるようにファイルを変換します.

program main
  implicit none
 
  integer(8)::nx,ny,nz,xx,yy,zz
  real(8),allocatable::x(:,:,:)
  real(8),allocatable::y(:,:,:)
  real(8),allocatable::z(:,:,:)
  real(8),allocatable::v(:,:,:)
 
  nx=50
  ny=50
  nz=50
 
  open(62,file="data.dat",status="old")!読み込みファイル
  open(64,file="x.txt",status="unknown")!出力ファイル
  open(65,file="y.txt",status="unknown")!出力ファイル
  open(66,file="z.txt",status="unknown")!出力ファイル
  open(67,file="v.txt",status="unknown")!出力ファイル

  allocate(x(-nx:nx,-ny:ny,-nz:nz))
  allocate(y(-nx:nx,-ny:ny,-nz:nz))
  allocate(z(-nx:nx,-ny:ny,-nz:nz))
  allocate(v(-nx:nx,-ny:ny,-nz:nz))
 
  do zz=-ny,ny
     do yy=-nx,nx
        do xx=-nz,nz
           read(62,*)x(xx,yy,zz),y(xx,yy,zz),z(xx,yy,zz),v(xx,yy,zz)
        end do
     end do
  end do
 
  do xx=-nx,nx
     write(64,*) x(xx,:,:)
     write(65,*) y(xx,:,:)
     write(66,*) z(xx,:,:)
     write(67,*) v(xx,:,:)
  end do
 
  deallocate(x)
  deallocate(y)
  deallocate(z)
  deallocate(v)
 
  close(62)
  close(64)
  close(65)
  close(66)
  close(67)
           
end program

これで, x.txt, y.txt, z.txt, v.txtが出力されます.
nx, ny, nzの値は, 範囲設定を変えた場合は適宜それに合わせて書き換えて下さい.
さて, 準備ができたので, MATLABを起動しましょう.

x.txt, y.txt, z.txt, v.txtの階層に移動して,
次のコマンドを入力してこれらのファイルを読み込みます.
次のコマンドはスクリプトファイルにしちゃうと楽ですね.

xn=50*2+1   %xのデータ数
yn=50*2+1   %yのデータ数
zn=50*2+1   %zのデータ数
x=load('x.txt')
x2=reshape(x,xn,yn,zn)
y=load('y.txt')
y2=reshape(y,xn,yn,zn)
z=load('z.txt')
z2=reshape(z,xn,yn,zn)
v=load('v.txt')
v2=reshape(v,xn,yn,zn)

コードの内容は,
まず, 先ほどのデータを読み込みます.
次に, 読み込んだデータをちょっと整理してます.

さて, あとは描くだけです！！
次のコマンドを入力していきましょう.

E=0.5 %Eが描きたい等値面の値
p = patch(isosurface(x2,y2,z2,v2,E));
isonormals(x2,y2,z2,v2,p)
set(p,'FaceColor','green','EdgeColor','none'); %ここで色変更可
daspect([1,1,1])
view(30,30); axis([-1.5 1.5 -1.5 1.5 -1.5 1.5]) %axisで表示範囲を設定
camlight ('headlight') %光の当てる方向
camlight ('right') %光の当てる方向
lighting gouraud

これで, 綺麗な楕円体が描けると思います.

hcu

こんな感じですかね.

もうちょっと頑張ったら, これの断面図をグロテスクに描いたりもできますが,
まぁアクセス数次第ですかね.

レッツ!!!等値面！！！！！！！

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

非線形Schrödinger方程式のソリトン解

コメントする

非線形シュレディンガー方程式
\(
\displaystyle i\frac{\partial \Psi}{\partial t}=-\frac{1}{2}\frac{\partial^2 \Psi}{\partial x^2}-g|\Psi|^2\Psi
\)
にはある解析解が存在します。それがソリトン(soliton)解と呼ばれるもので,上式のソリトン解は
\(
\displaystyle \Psi(x,t)=\sqrt{\Omega}\ {\rm sech}\left\{\sqrt{\Omega}\left(x\sqrt{g}-\frac{gV}{2}t\right)\right\}\cdot \exp\left\{i\frac{V\sqrt{g}}{2}x-i\frac{g}{2}\left(\frac{V^2}{4}-\Omega\right)t\right\}
\)
です。(\(g>0\),\(\Omega\):ソリトンの振幅、\(V\):ソリトンの速度に関するパラメータ、ソリトン自体の速度は\(V\sqrt{g}/2\))

[adsense1]

ソリトンの歴史的背景

「非線形」とは重ね合わせの原理が成り立たない系です。

1844年、スコットランドのJ.Scott-Russellによって孤立した波(solitary wave)を観測した事が報告されました J.Scott-Russellによる報告”Report on Waves”(リンク先のSR44.pdf, 16.3MB))。
当時の認識では、波は波動方程式で記述され、その波の速度\(v\)は\(v=f\lambda\)の元、一定である。だからパルス状の波は異なる波長の波の重ね合わせで書けているはずで、時間と共に分散していくはず。なのになぜ時間が経過しても孤立した波が存在できるのか?という事で大きな論争となりました。

60年後の1895年、オランダのKortewegとde Vriesによって”浅い水の波”を記述する非線形偏微分方程式(KdV方程式)が提出され、この方程式の特解として孤立波が与えられました。
孤立波は、

  1. 空間的に局在した波が、その性質(速さや形)を変えずに伝搬する
  2. 孤立波は互いの衝突に対して安定であり、各々の個別性を保持する

という性質を持つ非線形波動と定義されます[1]。
2番目の、粒子のような性質を持つことから、solitary に接頭語-on をつけ、soliton(ソリトン)と名づけられました。

その後、1981年に佐藤幹夫がソリトンの統一理論(佐藤理論やKP理論)を発表しました。
これによりソリトン方程式(ソリトンを記述し，かつ厳密に解ける方程式)に決着が付きました。
ソリトン方程式は非線形なのに厳密に解ける、可積分系である。

ソリトン方程式を解く方法は([4]を引用しますが)

上で指摘したように，logistic方程式が解けるからくりとソリトン方程式が解けるからくりはよく似ています．違いは，logistic方程式が変数変換一発で線形常微分方程式になってしまったのに対し，ソリトン方程式の場合は変数変換で双線形形式になり，双線形形式の解として行列式が現れ，行列式の中身に簡単な線形方程式が現れるというところです．しかし，離散化で保存するべき構造は明らかです．まず，解の中身の線形方程式を離散化し，行列式の構造をそのまま使って双線形形式を作る．最後に変数変換して非線形のレベルに戻ればよい．

となるそうです。

また、ソリトン方程式の特徴である、無限個の対称性(無限個の保存量)は、Gardner変換という変換をすることで証明できるそうです[5]。
これ以上はこの分野の専門家ではないので話せません。

ちなみに津波もソリトンの一つとみなせます。

ソリトン解が生まれるイメージ

なぜソリトン解が生まれるのでしょうか。
今、孤立した波(空間的に凸)を考えます。この時、

エネルギー的に安定になろうとして密度を均一にするために広がろうとする効果
と
粒子間を結び付ける引力相互作用(例えば水面だったら水と水との分子間力等)のため集まろうとする効果

のつり合いによって、丁度均衡が保たれるとき、このソリトン解が生まれます。

･･･実は、ソリトン解には2種類あります。
それは明るい(Bright)ソリトン解と暗い(Dark)ソリトン解です。
今まで話していたのは全て明るいソリトン解です。
暗いソリトン解とはどういったものでしょう。
暗いソリトン解とは、ある部分が空間的に凹んでいる、孤立した解です。

エネルギー的に安定になろうとして密度を均一にするためにその凹みを埋めようとする効果
と
粒子間の斥力相互作用のために粒子間を避けようとする効果

のつり合いによって、丁度均衡が保たれるとき、この暗いソリトン解が生まれます。
暗いソリトン解が生まれるのは斥力相互作用の時で、斥力相互作用を持つ系というのは、調べた限りでは量子力学のボーズ・アインシュタイン凝縮体で、暗いソリトンは渦ソリトンという形で現れるそうです。これ以上の具体例は分かりませんでした。もしも具体例を知っているという方は教えていただければ幸いです。
暗いソリトンの解析解は参考文献[1]の本に紹介されているので、それをご参考にしてください。

非線形シュレディンガー方程式におけるソリトン解

では本題の、非線形シュレディンガー方程式における(明るい)ソリトン解を考えましょう。
下の形の非線形シュレディンガー方程式を考えます。
\(
\displaystyle i\frac{\partial \Psi}{\partial t}=-\frac{1}{2}\frac{\partial^2 \Psi}{\partial x^2}-g|\Psi|^2\Psi
\)
ここで、\(\Psi=\Psi(x,t)\)で、\(g\)は正の値で相互作用の強さ(この場合、引力相互作用)を表します。

この非線形シュレディンガー方程式のソリトン解\(\Psi(x,t)\)は、
\(
\displaystyle \Psi(x,t)=\sqrt{\Omega}\ {\rm sech}\left\{\sqrt{\Omega}\left(x\sqrt{g}-\frac{gV}{2}t\right)\right\}\cdot \exp\left\{i\frac{V\sqrt{g}}{2}x-i\frac{g}{2}\left(\frac{V^2}{4}-\Omega\right)t\right\}
\)
であり、\(\Omega\)はソリトンの振幅の大きさ、\(V\)はソリトンの速度を決めるパラメータを表します。ソリトン自体の速度は\(V\sqrt{g}/2\)となります([2]を参考)。
また、\({\rm sech}(x)\)は双曲線関数の一種(双曲線正割と呼ばれる)であり、
\(
\displaystyle {\rm sech}(x)=\frac{1}{\rm cosh(x)}=\frac{2}{e^{x}+e^{-x}}
\)
を表します。

解析解のプロット

解析解をプロットします。gnuplotコードは下のほうに載せておきます。
\(g=2, V=1, \Omega=1\)とすると、以下の振る舞いが観測されます。
ここで紫はソリトン解の実部、緑は虚部、青は絶対値2乗を表します。
動画は1枚当たり原子単位系で0.1秒、合計で10秒間のシミュレーションです。
また、このソリトンの速度は\(V\sqrt{g}/2\sim 0.7071\)です。
soliton1

[adsense2]

gnuplotコード

非線形シュレディンガー方程式のソリトン解(解析解)を出力します。
gnuplot上で下のスクリプトを実行してください。
(ただし、gnuplot ver4.6以降に限ります。)

omega=1e0
V=1e0
g=2e0
x0=-5e0 # initial position

set xr[-10:10]
set yr[-1.5:1.5]
set samples 1000
set xl "x[a.u.]"

sech(x)=2e0/(exp(x)+exp(-x))
amp(x,t)=sqrt(omega)*sech(sqrt(omega)*((x-x0)*sqrt(g)-g*V*t*0.5e0))
phase(x,t)=V*sqrt(g)*x*0.5e0-g*t*0.5e0*(V*V*0.25e0-omega)
soliton(x,t)=amp(x,t)*exp({0e0,1e0}*phase(x,t))


#set term gif animate delay 10 optimize size 960,720
#set output 'movie.gif'
do for[i=0:100:2]{
   t=i*0.1e0
   plot abs(soliton(x,t))**2 lw 3 lc 1 lt 1 ti sprintf("|\psi|^2, t=%2.1f",t),\
        real(soliton(x,t)) lw 3 lc 2 lt 2 ti "Real",\
        imag(soliton(x,t)) lw 3 lc 3 lt 3 ti "imag"
}
#set out
#set terminal wxt enhanced

もっとソリトンについて知りたい方はまず参考文献[3]を読むことをお勧めします。
その後、[4]を読み、[5]を読み、[1]の本を読むのが良いと思われます。
[3],[4]は簡単な表現を用いてソリトンとその後の発展について記述されています。

参考文献

[1]和達三樹著 『非線形波動 (現代物理学叢書) 』岩波書店 (2000年) p.7
[2]和達三樹著 『非線形波動 (現代物理学叢書) 』岩波書店 (2000年) p.29

[3]ソリトンの数学 – Researchmap
[4]ソリトン ～ 不思議な波が運んできた，古くて新しい数学の物語 ～
[5]〔連載〕非線形波動―ソリトンを中心として―第５章　逆散乱法

↑画像のフォントはキユマヤ園様による数式フォント -びゅんびゅん→SSSです！

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

fortran90でファイルの存在を確かめる

コメントする

fortran90です。
ファイルが存在するかしないかを確かめるためのプログラムです。
組み込み関数”access”を用いると簡単に確認できます。

program main
  implicit none
  integer::access
 
  write(6,*)access( "./temp.d", " ")
 
  stop
end program main

上の文は、もしもファイル “./temp.d” が今いるディレクトリに存在するかどうかを返すものです。

また、access関数はファイルの有無(引数は空白文字” “)のみならず、そのファイルの読み取り権限(引数は”r”)、書き込み権限(引数は”w”)、実行権限(引数は”x”)を調べることができます。

ゼロが返ってきたらその権限(ファイルの有無)があるということです。
実行例は

 $ ls
main.f90  temp.d
 $ gfortran main.f90
 $ ./a.out
           0
 $

となります。intel®fortranコンパイラでももちろんできます。

[adsense1]

ここより以下のものはaccess関数を用いない方法です。

ファイルをopenできたら存在する、errorが返ってきたら存在しない、と判断します。
そのプログラムとサブルーチンはこんな感じ。

program main
  implicit none
  integer::ox
 
  call checkfile("temp.d",ox)
 
  write(6,*)ox
 
  stop
end program main

subroutine checkfile(fname,ox)
  implicit none
  character(*),intent(in)::fname
  integer,intent(out)::ox
 
  open(100,file=trim(fname),status='old',err=999)
  close(100)
  write(6,'(3A)')"file '",fname,"' exist"
  ox=1
  return
 
999 continue
  close(100)
  write(6,'(3A)')"file '",fname,"' don't exist"
  ox=0

  return
end subroutine checkfile

”ox”は〇×のつもりで使っています。
ファイル名”temp.d”をinputパラメータとして渡し、そのファイルが存在したらox=1を、存在しなかったらox=0を返します。

実行例は、

 $ ls
   temp.d main.f90
 $ gfortran main.f90
 $ ./a.out
   file 'temp.d' exist
           1
 $ gfortran main.f90
 $ ./a.out
   file 'temp1.d' don't exist
           0

となります。２回目の実行ではtemp.dをtemp1.dに変えて存在しない場合を表示しています。

[adsense2]

参考
ファイルの有無を調べる -知識の箱