sikino のすべての投稿

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

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

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と分かっていれば、筆算の計算は途中までは掛け算と引き算で行けます。そう考えるとまぁ多少は簡単なんですかね・・・？

参考文献

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

LaTeX

latexで数式全体を小さくする

2件のコメント

数式が長すぎてはみ出してしまう時の対処法

数式全体を小さくするには、数式を挟むように
\footnotesize
\normalsize
と書けばいいです。

\footnotesize
\begin{align}
g(x)=\sum_n c_n f_n(x)
\end{align}
\normalsize
実例

通常だと

\begin{align}
g(x)=\sum_n c_n f_n(x)
\end{align}

\footnotesize,\normalsizeで挟むと、

\footnotesize
\begin{align}
g(x)=\sum_n c_n f_n(x)
\end{align}
\normalsize

\tiny,\normalsizeで挟むと、

\tiny
\begin{align}
g(x)=\sum_n c_n f_n(x)
\end{align}
\normalsize

となります。

上から順に通常、footnotesize、tinyで挟んだ時を出力すると以下のようになります。

参考文献
Reducing font in equation -latex-community

gnuplot

gnuplotで点にラベルを付けて出力

コメントする

gnuplotで点とその点の値を付けたいとします。
数点だけならば、

set label "moji" at 4,0.2

とやれば点(4,0.2)に”moji”という文章が書かれます。

複数点の場合を考えます。
ここでは、データ”test.d”として、

1 0.1
2 0.2
3 0.3
4 0.4
5 0.5
6 0.6
7 0.7

を考えます。gnuplot単体でやることはできないと思います。
若干作業しなければなりません。

ファイル”test.plt”を生成し、中に以下の文を書きます。

unset label

!perl -ane 'print "set label sprintf(\"\%6.4f\",$F[1]) at $F[0]-0.05,$F[1]+0.04\n"' test.d > p.plt

load "p.plt"
plot "test.d" pt 7

とすると各点の位置からx軸方向に-0.05, y軸方向に+0.04ずれた位置にその点のy軸の値が出力されます。gnuplot上で上のファイルを読み込むと以下の出力を得ます。
point_label

参考文献
Basic statistics with gnuplot

gnuplot

gnuplotである1つの点のみを表示する

コメントする

gnuplotで複数行から成るデータの第n行目だけを点としてプロットしたい場合があります。

例えば、全100行から成るデータ(“data.d”)の第17行目だけをプロットしたいとします。
この時、

plot "data.d" every ::17:0:17:0 u 1:2

とすれば、第17行目だけが点としてプロットされます。

その他、例えばデータではなく、ある関数のf(x)のx=5.7の時だけの値をプロットしたい時は

plot sprintf("< echo ''%f %f''", 5.7,f(5.7))

とすればよいです。

例1

データ”test.d”として

1 0.1
2 0.2
3 0.3
4 0.4
5 0.5
6 0.6
7 0.7

を考えます。この時gnuplot上で、

plot "test.d" every ::5:0:5:0 pt 7

と行うと
pointevery
というグラフが得られます。

例2

sin(5.7)を考えます。gnuplot上で

plot sprintf("< echo ''%f %f''", 5.7,sin(5.7)) pt 7

とすれば
point_function
という画像を得ます。都合上同時にsin(x)も出力しています。

参考文献
　複数のデータファイルのある列の差分表示 -計算機は楽し
gnuplotのコマンドまとめ(ループとかeveryとか) -Linux, Mac, Emacsについての設定、覚え書き

弾道計算

上下幅が一定の時の軌道と、撃ち上げ撃ち下し時の軌道のずれについて

コメントする

今まではゼロインする距離を固定した時の比較を行ってきましたが、
通常軌道を見て適正ホップだ、とか決定するかと思います。
この条件ならば上下幅を一定にしなければ適切な比較とは言えないと思ったのでその場合の結果です。

また、山のフィールドで
斜面の下から上に向かって打つとき、
斜面の上から下へ向かって打つとき
の着弾点は上にずれることを示します。

上下幅が一定の時の軌道

弾道の上下の幅を、0.25gのBB弾で0.9Jで射出した時の弾道の上下幅を基準にして重さを変えた時、弾道がどうなるかを計算すると以下の通りになります。
変更するBB弾の重さは販売されているものを基準とし、
0.12g, 0.20g, 0.21g, 0.23g, 0.25g, 0.28g, 0.30g
にしています。
グラフ中の点は0.05秒ごとの位置をプロットしており、図中の上の数本の線は水平距離のみを見た時のものです。
orbit_updown

静止画では
重さと軌道の関係_c
のようになります。

飛距離に対して到達速度は0.12g以外ではほとんど変わりません。
30m付近でおおよそ同じ時間に到達していることが分かります。
それ以降は重いほど早く到達し、かつ遠くまで飛んでいることが分かるかと思います。

やはり銃が許す限り、BB弾は重ければ重いほど良い、ということです。

撃ち上げ、撃ち下し時の軌道のずれについて

山のフィールドにおいて、斜面下から上、または上から下に向かって撃つ時、軌道はどう変わるのでしょうか。
答えは、照準よりも上に飛んでいきます。

以下の画像の状況を考えます。
角度_c

エネルギーは0.90J, 重さは0.25gを考えます。
下から上に向かって打つときの角度を正にとって\(+\theta\)
上から下に向かって打つときの角度を負にとって\(-\theta\)と表すことにします。
また、平坦な場所でゼロイン50mで最小の振れ幅になる軌道で比較します。

それぞれの角度\(\theta\)の向きで射出した後、撃った角度\(\theta\)で座標を回転させて主観的な視点で見た時にどうなるかを考えます。

すると、以下の結果を得ます。
撃ち上げ撃ち下し時_軌道のずれ_c
横軸は水平の距離ではなく、その角度で見た時の相対的な距離です。
\(\theta=0\)の時とそれぞれの角度の時を比較すると、角度をつける場合はほとんどの場合で上向きにずれることが分かります。
すなわち、角度がある時はシューティングレンジで合わせた照準よりも上に着弾するのです。

何故でしょうか。ちゃんと理由があります。
上向きにでも下向きにでも、角度があって撃つとき、BB弾に掛かる見かけ上の重力が減少します。
gravity_updown_c
なので角度がある時は重力が小さく感じるので、角度が無いところで調節した時よりも上の方に飛んでいくことになります。
なので、山のフィールドで斜面があるときは、若干下方向を、すなわち相手の膝や腰辺りを狙えばいいのです。

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

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で現在開いてるバッファを再読込する