ファイル読み込み

ファイルの読み込みに関するサブルーチンを記述します。

 do i=1,4
     do j=1,83
        write(30,*)i,j,i*j
     end do
     write(30,*)
  end do

によってファイル”fort.30″が作られたとします。
今、fort.30を読み込んで
“i”の数”4”

“j”の数”83″
を取り出したいとします。

read文を使って読み込みますが、そのまま読み込むと空白部分を読み込んだり、読み込まないだったりします。
ここでは上のdoループによって作られたファイルの場合に使えるものを想定します。

ここでは、
大きな塊を表す数である”4″をblock,
塊の中の要素数を表す”83″をelement
と呼ぶことにします。

この問題を解く戦略は、ファイルを読み込む2種類の方法です。
1つは本当の行数(空行含む)を数え、もう1つは空行を飛ばして読み込む方法です。
本当の行数を数えるサブルーチンは下の”linecount”
であり、
空行を読み飛ばして行数を数えるサブルーチンは下の”linecount_eff”
です。
この二つと、一番下の行に追加される余分な1行を対処するために”breaklinecheck”というルーチンを使います。

これによってblockとelementを出力するサブルーチン”blockelement”を作っています。

下のプログラムを実行すると

$ gfortran main.f90
$ ./a.out
 ===Apply for fortran file will done===
 Nblock     ::         4
 Nelement   ::        83

という結果が得られるかと思います。

program main
  !developer => sikinote
  !date      => 2015/03/31
  implicit none
  integer::NBlock,Nelement
  character(48)::filename
 
  filename='./fort.30'
  call blockelement(filename,Nblock,Nelement)

  write(6,'(A,i10)')" Nblock     ::",Nblock
  write(6,'(A,i10)')" Nelement   ::",Nelement

  stop
end program
!===================================
   
subroutine blockelement(filename,Nblock,Nelement)
  !developer => sikino
  !date      => 2015/03/31
  implicit none
  character(*),intent(in)::filename
  integer,intent(out)::Nblock,Nelement
  integer::c1,c2
 
  call linecount(c1,filename)
  call linecount_eff(c2,filename)
  call breaklinecheck(c1,c2)
  Nblock=c1-c2+1
  Nelement=c2/Nblock
 
  return
end subroutine blockelement
!------------------------------
subroutine linecount(c,filename)
  implicit none
  integer,intent(out)::c
  character(*),intent(in)::filename

  integer::ier
  character(len_trim(filename))::fname
 
  fname=trim(filename)
  c=0
  open(100,file=fname,status='old',iostat=ier,err=990)
  do while(.true.)
     read(100,*,end=999)
     c=c+1
  enddo
999 continue
  close(100)
   
  return
!+-----------------------+
990 write(6,'(A)')"!!!!==error when open file",trim(fname),"info-->",ier
  write(6,*)"======program stop at linecount"
  stop

end subroutine linecount
!--------------------------------
subroutine linecount_eff(c,filename)
  implicit none
  integer,intent(out)::c
  character(*),intent(in)::filename

  integer::ier
  character(100)::cc
  character(len_trim(filename))::fname
 
  fname=trim(filename)

  c=0
  open(100,file=fname,status='old',iostat=ier,err=990)
  do while(.true.)
     read(100,*,end=998)cc
     if(len_trim(cc).gt.0)c=c+1
  enddo

998 continue
  close(100)
  return

990 write(6,'(A)')"!!!!==error when open file",trim(fname),"info==>",ier
  write(6,*)"======program stop at linecount_eff"
  stop

end subroutine linecount_eff
!-------------------------------------
subroutine breaklinecheck(c1,c2)
  implicit none
  integer,intent(inout)::c1
  integer,intent(in)::c2
  integer::Nb
 
  Nb=c1-c2+1
  if(Nb.eq.2.or.mod(c2,Nb).ne.0)then
     write(6,*)"===Apply for fortran file will done==="
     c1=c1-1
     Nb=c1-c2+1
     if(mod(c2,Nb).ne.0)then
        write(6,*)"line is different(may be last break)"
        write(6,*)"program stop at subroutine __breaklinecheck__"
        stop
     end if
  endif
 
  return
end subroutine breaklinecheck

データを読み込むには?


さて、ブロックの数と要素の数が上のサブルーチンを使うことにより求められることがわかりました。
実際にデータを配列に代入するためにはどうすればいいんでしょう?
型に応じて使うサブルーチンを変えます。
その手続きは下のモジュールを記述することでokです。これを書いた上で、
メインプログラムを以下のように書きます。そうすれば配列x(,)とy(,)に値がはいります。

program main
  use read1
  implicit none
  integer::NBlock,Nelement,i,j
  character(48)::filename
 
  double precision,allocatable::x(:,:),y(:,:)
 
  filename='./fort.30'
  call blockelement(filename,Nblock,Nelement)

  write(6,'(A,i10)')" Nblock     ::",Nblock
  write(6,'(A,i10)')" Nelement   ::",Nelement
 
  allocate(x(1:Nblock,1:Nelement),y(1:Nblock,1:Nelement))
  call read_filedata(size(y,1),size(y,2),x,y,filename)
 
  do i=1,Nblock
     do j=1,Nelement
        write(11,*)x(i,j),y(i,j)
     enddo
     write(11,*)
  enddo
 
  stop
end program

総称名を用いる場合の手続き(read_filedata())

module read1
  implicit none
  interface read_filedata
     module procedure &
          ! dx1 -> double precision array, x(:)
          ! cy2 -> complex array, y(:,:)
          ! xyy -> coloum of file, | x y y |
          read_dx0_dy1_xy, &
          read_dx0_cy1_xyy, &
          read_dx1_dy1_xy, &
          read_dx1_cy1_xyy, &
          read_dx1_dy2_xy, &
          read_dx1_cy2_xyy, &
          read_dx2_dy2_xy, &
          read_dx2_cy2_xyy
  end interface read_filedata
contains
  subroutine read_dx0_dy1_xy(Ne,y,place,col12)
    integer,intent(in)::Ne,col12
    character(*),intent(in)::place
    double precision,intent(out)::y(1:Ne)

    integer::i,ier
    double precision::a,b
    character(len_trim(place))::fn

    y=0d0
    fn=trim(place)

    open(28,file=fn,status='old',iostat=ier,err=977)
    if(col12.eq.1)then
       do i=1,Ne
          read(28,*)a,b
          y(i)=a
       enddo
    elseif(col12.eq.2)then
       do i=1,Ne
          read(28,*)a,b
          y(i)=b
       enddo
    else
       go to 977
    endif
             
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx0_dy1_11"
    stop
  end subroutine read_dx0_dy1_xy

  subroutine read_dx0_cy1_xyy(Ne,y,place)
    integer,intent(in)::Ne
    character(*),intent(in)::place
    complex(kind(0d0)),intent(out)::y(1:Ne)

    integer::i,ier
    double precision::a,b,c
    character(len_trim(place))::fn

    y=0d0
    fn=trim(place)

    open(28,file=fn,status='old',iostat=ier,err=977)
    do i=1,Ne
       read(28,*)a,b,c
       y(i)=dcmplx(b,c)
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx0_cy1_xyy"
    stop
  end subroutine read_dx0_cy1_xyy

  subroutine read_dx1_dy1_xy(Ne,x,y,place)
    integer,intent(in)::Ne
    character(*),intent(in)::place
    double precision,intent(out)::x(1:Ne),y(1:Ne)

    character(len_trim(place))::fn
    double precision::a,b
    integer::i,ier

    x=0d0; y=0d0
   
    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Ne
       read(28,*)a,b
       x(i)=a
       y(i)=b
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx1_dy1_xy"
    stop
  end subroutine read_dx1_dy1_xy
 
  subroutine read_dx1_cy1_xyy(Ne,x,y,place)
    integer,intent(in)::Ne
    character(*),intent(in)::place
    double precision,intent(out)::x(1:Ne)
    complex(kind(0d0)),intent(out)::y(1:Ne)

    integer::i,ier
    double precision::a,b,c
    character(len_trim(place))::fn

    y=0d0
    fn=trim(place)

    open(28,file=fn,status='old',iostat=ier,err=977)
    do i=1,Ne
       read(28,*)a,b,c
       x(i)=a
       y(i)=dcmplx(b,c)
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx1_cy1_xyy"
    stop
  end subroutine read_dx1_cy1_xyy

  subroutine read_dx1_dy2_xy(Nb,Ne,x,y,place)
    integer,intent(in)::Nb,Ne
    double precision,intent(out)::x(1:Ne),y(1:Nb,1:Ne)
    character(*),intent(in)::place
   
    integer::i,j,ier
    character(len_trim(place))::fn
    double precision::a,b
   
    x=0d0; y=0d0
   
    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Nb
       do j=1,Ne
          read(28,*)a,b
          if(i.eq.1)x(j)=a
          y(i,j)=b
       enddo
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx1_dy2_xy"
    stop
  end subroutine read_dx1_dy2_xy
 
  subroutine read_dx1_cy2_xyy(Nb,Ne,x,y,place)
    integer,intent(in)::Nb,Ne
    double precision,intent(out)::x(1:Ne)
    complex(kind(0d0)),intent(out)::y(1:Nb,1:Ne)
    character(*),intent(in)::place

    integer::i,j,ier
    character(len_trim(place))::fn
    double precision::a,b,c

    x=0d0; y=0d0

    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Nb
       do j=1,Ne
          read(28,*)a,b,c
          if(i.eq.1)x(j)=a
          y(i,j)=dcmplx(b,c)
       enddo
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx1_cy2_12"
    stop
  end subroutine read_dx1_cy2_xyy

  subroutine read_dx2_dy2_xy(Nb,Ne,x,y,place)
    integer,intent(in)::Nb,Ne
    double precision,intent(out)::x(1:Nb,1:Ne),y(1:Nb,1:Ne)
    character(*),intent(in)::place
   
    integer::i,j,ier
    character(len_trim(place))::fn
    double precision::a,b
   
    x=0d0; y=0d0
   
    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Nb
       do j=1,Ne
          read(28,*)a,b
          x(i,j)=a
          y(i,j)=b
       enddo
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx2_dy2_11"
    stop
  end subroutine read_dx2_dy2_xy

  subroutine read_dx2_cy2_xyy(Nb,Ne,x,y,place)
    integer,intent(in)::Nb,Ne
    double precision,intent(out)::x(1:Nb,1:Ne)
    complex(kind(0d0))::y(1:Nb,1:Ne)
    character(*),intent(in)::place
   
    integer::i,j,ier
    character(len_trim(place))::fn
    double precision::a,b,c
   
    x=0d0; y=dcmplx(0d0,0d0)
   
    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Nb
       do j=1,Ne
          read(28,*)a,b,c
          x(i,j)=a
          y(i,j)=dcmplx(b,c)
       enddo
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx2_cy2_12"
    stop
  end subroutine read_dx2_cy2_xyy
end module read1

ラグランジュの未定乗数法

忘れやすいラグランジュの未定乗数法のメモです。
一通り学んだ人が使い方を思い出す、という状況を想定しています。

変数が独立な場合


3変数x,y,zが独立(xが変化してもy,zは変化しない、\vec{x}\cdot\vec{y}=0 (y,zも同様))で、その関数f(x,y,z)の極値は
\displaystyle \frac{\partial f}{\partial x}=0,\ \ \frac{\partial f}{\partial y}=0,\ \ \frac{\partial f}{\partial z}=0 \ \ \cdots (1)
を連立させて解くことで得られます。

変数が独立ではない(従属な)場合


変数x,y,zの間に関係式
g(x,y,z)=c, \ \mbox{$c$は定数} \ \ \cdots (2)
という条件がある場合、

  1. 式(2)が(例えば)zについて解けるならば、f(x,y,z(x,y))として2変数x,yの極値問題として解ける。
  2. z=z(x,y)の形に書けない場合 →ラグランジュの未定乗数法を使う
  1.  の場合、”関数f(x,y,z)が極値をとる”、とは
    \displaystyle df=\frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial y}dy+\frac{\partial f}{\partial z}dz=0 \ \ \cdots (3)
    である点(極値点)である。
    →極値点からx,y,zを任意の微小量dx,dy,dzだけ変化させても関数f(x,y,z)の変化分であるdfは1次の範囲でゼロです。
    故に任意のdx,dy,dzについて(3)が成立する、よって(1)が導かれることになります。
  2.  の場合、dx,dy,dz独立に選ぶことができません。その取り方は条件(2)に従います。
    極値点まわりで(2)が成立しているならば、
    g(x+dx,y+dy,z+dz)=g(x,y,z)\ \ \cdots (4)
    を満たすようなdx,dy,dzの変化しか許されないことになります。この条件からdx,dy,dzの動かし方は、
    \displaystyle \frac{\partial g}{\partial x}dx+\frac{\partial g}{\partial y}dy+\frac{\partial g}{\partial z}dz=0 \ \ \cdots (5)
    に制限されます。
    条件(5)をdzについて変形すると、
    \displaystyle dz=-\frac{\frac{\partial g}{\partial x}dx+\frac{\partial g}{\partial y}dy}{\frac{\partial g}{\partial z}} \ \ \cdots (6)
    となります。式(3)へ式(6)を代入すると

    \begin{align} \displaystyle df &=\left( \frac{\partial f}{\partial x}-\frac{\frac{\partial g}{\partial x}}{\frac{\partial g}{\partial z}}\frac{\partial f}{\partial z}\right)dx +\left( \frac{\partial f}{\partial y}-\frac{\frac{\partial g}{\partial y}}{\frac{\partial g}{\partial z}}\frac{\partial f}{\partial z}\right)dy=0 \\ &=\left( \frac{\partial f}{\partial x}-\lambda\frac{\partial g}{\partial x}\right)dx +\left( \frac{\partial f}{\partial y}-\lambda\frac{\partial g}{\partial y}\right)dy=0 \ \ \cdots (7) \end{align}

    ここで
    \displaystyle \lambda=\frac{\frac{\partial f}{\partial z}}{\frac{\partial g}{\partial z}}\ \ \cdots (8)
    と置きました。この\lambdaはラグランジュの未定乗数と呼ばれます。
    未定乗数という所以は、この\lambdaをあらわに決める必要はなく、決まらない定数のままでも極値点を求めることができるという事を表しています。
    今、dx,dyは独立であるので(※1)、その係数は0になるはずです。故に、式(7)と式(8)より、

    \begin{eqnarray} \left\{ \begin{aligned} \frac{\partial f}{\partial x}-\lambda\frac{\partial g}{\partial x}&=0 \\ \frac{\partial f}{\partial y}-\lambda\frac{\partial g}{\partial y}&=0 \\ \frac{\partial f}{\partial z}-\lambda\frac{\partial g}{\partial z}&=0 \end{aligned} \right. \end{eqnarray}

    が導けます。変形をすれば、条件式(2)も含めて、

    \begin{eqnarray} \left\{ \begin{aligned} \frac{\partial}{\partial x}\left(f-\lambda g\right)&=0 \\ \frac{\partial}{\partial y}\left(f-\lambda g\right)&=0 \\ \frac{\partial}{\partial z}\left(f-\lambda g\right)&=0 \\ g(x,y,z)&=c \end{aligned} \right. \end{eqnarray} \ \ \ \cdots (9)

    と書くことができます。未知の変数はx,y,z,\lambdaの4つで、方程式は4本なので解くことができます。
    この式が言うことは、束縛条件g(x,y,z)=cがあった場合、関数f-\lambda gを考えて、その極値を求めればよいことを表しています。

[adsense1]

まとめ


条件g(x,y,z)=cがある関数f(x,y,z)の極値問題は、
\tilde{f}=f-\lambda g \ \ \cdots (10)
という新たな関数\tilde{f}を考えるとx,y,zが独立に変化するものと考えて\tilde{f}(x,y,z)の極値問題を考えればよい、となります。

例題 ~楕円に内接する長方形の面積の最大値を求める~


楕円の方程式は
\displaystyle \frac{x^2}{a^2}+\frac{y^2}{b^2}=1
であり、この方程式の許すx,yを満たしながら、内接する長方形の面積
S(x,y)=4xy
を最大にするx,yを求めます。

楕円面積

これは、関数f(x,y)=4xyの極値を\displaystyle g(x,y)=\frac{x^2}{a^2}+\frac{y^2}{b^2}=1という条件下で解く、という意味になります。
まず、\tilde{S}を式(10)より求めます。
\begin{align} \tilde{S}&=S-\lambda g \\ &=4xy-\lambda\left(\frac{x^2}{a^2}+\frac{y^2}{b^2}\right) \end{align}
式(9)より、
\begin{eqnarray} \left\{ \begin{aligned} \frac{\partial \tilde{S}}{\partial x}&=0 \\ \frac{\partial \tilde{S}}{\partial y}&=0 \end{aligned} \right. \end{eqnarray}
を考えればいいので、
\begin{eqnarray} \left\{ \begin{aligned} \frac{\partial \tilde{S}}{\partial x}=4y-\lambda \frac{2x}{a^2}&=0 \ \ \cdots (i)\\ \frac{\partial \tilde{S}}{\partial y}=4x-\lambda \frac{2y}{b^2}&=0 \ \ \cdots (ii)\\ \frac{x^2}{a^2}+\frac{y^2}{b^2}&=1 \ \ \cdots (iii) \end{aligned} \right. \end{eqnarray}
を満たす未知の変数x,y,\lambdaSの極値となっています。
(i)と(ii)より\lambdaを消去すると
\begin{align} (i) &\rightarrow \lambda=\frac{2y}{x}a^2 \\ (ii) &\rightarrow 4x-\left(\frac{2y}{x}a^2\right)\frac{2y}{b^2}=0 \end{align}
なので
\displaystyle \frac{x^2}{a^2}=\frac{y^2}{b^2}
(iii)に代入して
\displaystyle 2\frac{x^2}{a^2}=1
より
\displaystyle x=\pm\frac{a}{\sqrt{2}},\ \ y=\pm\frac{b}{\sqrt{2}}
のときS(x,y)が極値を取ることがわかります。

ちなみに、この時の長方形の面積S_{max}S_{max}=2abであり、
これは楕円の面積S=\pi ab\frac{2ab}{\pi ab}\sim 0.64となり、約64%を占めていることになります。

[adsense2]

※1
x,yが独立であるのは、x,y,zをつなぐ1本の条件式g(x,y,z)=cによって消え得る変数は1つだけであるためです。

参考


小野寺 嘉孝著 『物理のための応用数学』裳華房(1988)p.6~10

振り子

振り子の角度が小さいとき、また大きいときの振り子の取扱いです。
ラグランジアンLとハミルトニアンHから出発して運動方程式を求める方法を記述します。
座標は図のようにとります。
振り子のラグランジアンとハミルトニアン2
式で表せば、
\begin{eqnarray} \left\{ \begin{aligned} x&=r \sin{\theta} \\ y&=-r \cos{\theta} \end{aligned} \right. \end{eqnarray}
であり、1階微分は
\begin{eqnarray} \left\{ \begin{aligned} \dot{x}&=\dot{r} \sin{\theta}+r\dot{\theta}\cos{\theta} \\ \dot{y}&=\dot{r} \cos{\theta} -r\dot{\theta}\sin{\theta} \end{aligned} \right. \end{eqnarray}
です。

[adsense1]


ラグランジアンLから出発


上の座標の元、デカルト座標(x,y)での振り子のラグランジアンL(x,\dot{x},y,\dot{y},t)
運動エネルギーK=\frac{1}{2}m(\dot{x}^2+\dot{y}^2), ポテンシャルエネルギーV=V(x,\dot{x},y,\dot{y},t) より、
\begin{align} L(x,\dot{x},y,\dot{y},t)&=K-V \\ &=\frac{1}{2}m(\dot{x}^2+\dot{y}^2)-V(x,\dot{x},y,\dot{y},t) \end{align}
と書けます。新たな座標系でのラグランジアンL(r,\dot{r},\theta,\dot{\theta},t)\begin{eqnarray} \left\{ \begin{aligned} x&=r \sin{\theta} \\ y&=-r \cos{\theta} \end{aligned} \right. \end{eqnarray}に従い座標変換をし、
\displaystyle L(r,\dot{r},\theta,\dot{\theta},t)=\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2)-V(r,\dot{r},\theta,\dot{\theta},t)
とすることで新しい座標でのラグランジアンが得られます。
ここまでは一般的な話です。

今、もしも振り子をつないでいる棒が時間変化しない、言い換えれば束縛条件x^2+y^2=l^2,(ここでlは定数)の元で考える場合、ポテンシャルエネルギーVは
\begin{align} V&=V(x,\dot{x},y,\dot{y},t)=-mgy \\ &=V(r,\dot{r},\theta,\dot{\theta},t)=mgr\cos\theta \end{align}
であり、\dot{l}=0なので
\displaystyle L(\theta,\dot{\theta},t)=\frac{1}{2}ml^2\dot{\theta}^2-mgl\cos\theta
となります。

運動方程式を求めるにはラグランジュの運動方程式
\displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}}\right)-\frac{\partial L}{\partial \theta}=0
に代入し、計算します。
\begin{align} \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}}\right)-\frac{\partial L}{\partial \theta}&=0 \\ \frac{d}{dt}(ml^2\dot{\theta})+mgl\sin{\theta}&=0 \\ ml^2\ddot{\theta}&=-mgl\sin\theta \\ \ddot{\theta}&=-\frac{g}{l}\sin\theta \end{align}
となり\thetaに関する微分方程式が導けました。

ハミルトニアンHから出発


ハミルトニアンHはラグランジアンLと一般化座標q_iと一般化運動量p_iを用いて、
\displaystyle H(\vec{p},\vec{q},t)\equiv\sum_i p_i\dot{q}_i-L(\vec{q},\vec{\dot{q}},t)
と定義されます。また、一般化運動量p_iは一般化座標q_iを用いて
\displaystyle p_i=\frac{\partial L(q_i,\dot{q}_i,t)}{\partial \dot{q_i}}
という関係で定義されます。


デカルト座標での一般化運動量p_xp_y
\begin{eqnarray} \left\{ \begin{aligned} p_x&=\frac{\partial L(x,\dot{x},y,\dot{y},t)}{\partial \dot{x}}=m\dot{x} \\ p_y&=\frac{\partial L(x,\dot{x},y,\dot{y},t)}{\partial \dot{y}}=m\dot{y} \end{aligned} \right. \end{eqnarray}
となるため、デカルト座標でのハミルトニアン
\begin{align} H(x,p_x,y,p_y,t)&=p_x\dot{x}+p_y\dot{y}-\left\{\frac{1}{2}m(\dot{x}^2+\dot{y}^2)-V(x,\dot{x},y,\dot{y},t)\right\} \\ &=\frac{p_x^2}{m}+\frac{p_y^2}{m}-\left\{\frac{1}{2}m(\frac{p_x^2}{m^2}+\frac{p_y^2}{m^2})-V(x,\dot{x},y,\dot{y},t)\right\} \\ &=\frac{1}{2m}(p_x^2+p_y^2)+V(x,\dot{x},y,\dot{y},t) \end{align}
と求められます。

新たな座標系でのは一般化運動量p_r, p_{\theta}
\begin{eqnarray} \left\{ \begin{aligned} p_r&=\frac{\partial L(r,\dot{r},\theta,\dot{\theta},t)}{\partial \dot{r}}=m\dot{r} \\ p_{\theta}&=\frac{\partial L(r,\dot{r},\theta,\dot{\theta},t)}{\partial \dot{\theta}}=mr^2\dot{\theta} \end{aligned} \right. \end{eqnarray}
なので、新たな座標系でのハミルトニアンH(r,p_r,\theta,p_{\theta},t)
\begin{align} H(r,p_r,\theta,p_{\theta},t)&=p_r\dot{r}+p_{\theta}\dot{\theta}-\left\{\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2)-V(r,\dot{r},\theta,\dot{\theta},t) \right\} \\ &=\frac{p_r^2}{m}+\frac{p_{\theta}^2}{mr^2}-\left\{ \frac{p_r^2}{2m}+r^2\frac{p_{\theta}^2}{2mr^4} – V(r,\dot{r},\theta,\dot{\theta},t) \right\} \\ &=\frac{1}{2m}\left(p_r^2+\frac{p_{\theta}^2}{r^2}\right) + V(r,\dot{r},\theta,\dot{\theta},t) \end{align}
と求められます。

今、もしも振り子をつないでいる棒が時間変化しない、言い換えれば束縛条件x^2+y^2=l^2, (ここでlは定数) の元で考える場合r=lを代入し、ポテンシャルエネルギーV
\begin{align} V&=V(x,\dot{x},y,\dot{y},t)=-mgy \\ &=V(l,\dot{l},\theta,\dot{\theta},t)=mgl\cos\theta \end{align}
と求められます。今、p_l=m\dot{l}=0なのでハミルトニアンは
\displaystyle H(\theta,p_{\theta},t)=\frac{p_{\theta}^2}{2ml^2}+mgl\cos{\theta}
となります。
ハミルトンの正準運動方程式
\begin{eqnarray} \left\{ \begin{aligned} \frac{\partial H}{\partial q_i}&=-\dot{p}_i \\ \frac{\partial H}{\partial p_i}&=\dot{q}_i \end{aligned} \right. \end{eqnarray}
に代入すれば
\begin{eqnarray} \left\{ \begin{aligned} \frac{d p_{\theta}}{dt}&=-mgl\sin\theta \\ \frac{d \theta}{dt}&=\frac{p_{\theta}}{ml^2} \end{aligned} \right. \end{eqnarray}
より、p_{\theta}を消去すれば、
運動方程式
\displaystyle \frac{d^2\theta}{dt^2}=-\frac{g}{l}\sin\theta
と導けます。

振り子の運動方程式


今、解くべき方程式は
\displaystyle \frac{d^2\theta}{dt^2}=-\frac{g}{l}\sin\theta
になります。

振り子の角度が小さいとき


まず、よく近似される角度が小さい時だけを議論するとします。
\sin\thetaはテイラー展開より、
\displaystyle \sin{\theta}= \theta-\frac{\theta^3}{3!}+\frac{\theta^5}{5!}-\frac{\theta^7}{7!}+\cdots =\sum_{n=1}^{\infty}(-1)^{n+1}\frac{\theta^{2n-1}}{(2n-1)!}
であるため、\theta\ll 1である場合は初めの1項だけ考慮し、
\sin{\theta}\sim \theta
と近似します。すると、
\displaystyle \frac{d^2\theta}{dt^2}=-\frac{g}{l}\theta
となります。

通常は解\theta(t)C_1\sin(\sqrt{\frac{g}{l}}t)+C_2\cos(\sqrt{\frac{g}{l}}t)とか置いてしまいますが、今回は数学的な手順に従って導出します。

まず、\displaystyle\frac{d}{dt}\left\{\left(\frac{d\theta}{dt}\right)^2\right\}という量を計算すると、
\displaystyle \frac{d}{dt}\left\{\left(\frac{d\theta}{dt}\right)^2\right\}=2\frac{d\theta}{dt}\cdot \frac{d^2\theta}{dt^2}
という関係があることが分かります。これを利用すると2階微分を1階微分に落とすことができます。
まず、運動方程式の両辺に\displaystyle 2\frac{d\theta}{dt}を掛け、整理していきます。すなわち、
\begin{align} 2\frac{d \theta}{dt}\cdot\frac{d^2 \theta}{dt^2}=-\frac{g}{l}2\frac{d \theta}{dt}\cdot\theta \\ \frac{d}{dt}\left\{\left(\frac{d\theta}{dt}\right)^2\right\}=\frac{d}{dt}\left\{-\frac{g}{l}\theta^2\right\} \\ \left(\frac{d\theta}{dt}\right)^2=-\frac{g}{l}\theta^2+C_0 \end{align}
となります。ここで\left(\frac{d\theta}{dt}\right)^2\thetaが虚数をとらない限り必ず正であることから、
-\pi\lt \theta\lt \piの範囲で(ただし、\sin\theta\sim \thetaの成り立つ近似の領域が\thetaが小さい場合だけなので、まぁまぁ可能な範囲、として考えます。)
どんなに変化しても積分定数C_0はそれより大きい値(\left(\frac{d\theta}{dt}\right)^2\geq 0にする値、C_0は必ず正である。)でなければなりません。
それを考えながら左辺の2乗を取り払い、
\begin{align} \frac{d\theta}{dt} = \pm \sqrt{C_0-\frac{g}{l}\theta^2} \\ \pm \int \frac{1}{\sqrt{C_0-\frac{g}{l}\theta^2}}d\theta=\int dt \end{align}
となります。左辺を計算するため変数変換
\displaystyle \sqrt{\frac{g}{l}}\theta=\sqrt{C_0}\sin\phi
を行います。
\displaystyle \frac{d \theta}{d\phi}=\sqrt{\frac{l}{g}C_0}\cos\phi
なので
\begin{align} \displaystyle \int \frac{1}{\sqrt{C_0-\frac{g}{l}\theta^2}}d\theta &=\displaystyle \int \frac{1}{\sqrt{C_0}\sqrt{1-\sin^2{\phi}}}\cdot\sqrt{\frac{l}{g}C_0}\cos{\phi}d\phi \\ &=\sqrt{\frac{l}{g}}\phi \\ &=\sqrt{\frac{l}{g}}\sin^{-1}\left(\sqrt{\frac{g}{lC_0}}\theta\right) \end{align}
と計算できるので、
\begin{align} \sqrt{\frac{l}{g}}\sin^{-1}\left(\sqrt{\frac{g}{lC_0}}\theta\right)=t+C_1 \\ \theta(t)=\sqrt{\frac{l}{g}C_0}\sin\left(\pm \sqrt{\frac{g}{l}}(t+C_1)\right) \end{align}
と書けます。C_0,C_1は単なる積分定数であるためにそれ自体の意味はなく、また、\sinの中の±はC_1を調節すれば関係なくなるので、それらを省略すれば
\displaystyle \theta(t)=C_0\sin\left(\sqrt{\frac{g}{l}}t+C_1\right)
となります。角周波数は
\displaystyle \omega=\sqrt{\frac{g}{l}}
となります。

角度が小さいときはこれでお仕舞いにします。

振り子の角度が大きいとき


楕円関数というものが計算を行っていくと出てきます。そこに至るまでの手順は振り子の角度が小さいときの手順に途中まで同じです。
対象になる問題は
\displaystyle \frac{d^2\theta}{dt^2}=-\frac{g}{l}\sin\theta
です。ここで\sqrt{\frac{g}{l}}=\omegaとおいて、
\displaystyle \frac{d^2\theta}{dt^2}=-\omega^2\sin\theta
と書き直します。
振り子の角度が小さいときと同じようにまず2階微分を1階微分に焼直すため、両辺に\displaystyle 2\frac{d\theta}{dt}を掛けて、
\begin{align} 2\frac{d \theta}{dt}\cdot\frac{d^2 \theta}{dt^2}=-\omega^2\frac{d \theta}{dt}\cdot \sin\theta \\ \frac{d}{dt}\left\{\left(\frac{d\theta}{dt}\right)^2\right\}=\frac{d}{dt}\left\{2\omega^2\cos\theta\right\} \\ \left(\frac{d\theta}{dt}\right)^2=2\omega^2\cos(\theta)+C_0 \end{align}
ここで左辺はゼロ以下になることはないため、そうなるように積分定数C_0を決めます。
ここでの積分定数C_0の持つ物理的な意味をちょっと考えます。
振り子には大きく分けて2種類の運動が存在します。1つは\thetaが小さいときの運動の延長線上にある、最下点を中心として振動する運動であり、もう1つはぐるんぐるん回る運動です。
ちょっと違う言い方をすれば、1つ目の運動は必ず粒子の速度の大きさ|v|=l|\frac{d \theta}{dt}|=0になる時間が存在する運動を表し、2つ目の運動は常に速度の大きさ|v|がゼロではない有限の値を持つ運動を表します。
この運動の区別がこの積分定数C_0に依存するのです。

何はともあれ、数学的な話をしていけば2つの運動を記述する微分方程式はどちらも\displaystyle \frac{d^2\theta}{dt^2}=-\omega^2\sin\thetaであるので、これが解ければ自然と両方の運動を記述する答えが得られるはずです。その区別は後で考えましょう。

つまり、言いたかったことはC_0はいつも2\omega^2\cos(\theta)以上ですよ、ということを言いたかったんです。
よって、左辺のルートを取り払います。
\begin{align} \left(\frac{d\theta}{dt}\right)^2 &=2\omega^2\cos(\theta)+C_0 \\ \frac{d\theta}{dt} &=\pm\sqrt{C_0+2\omega^2\cos(\theta)} \end{align}
天下り的ですが、ここで
\displaystyle \cos\theta=1-2\sin^2{\frac{\theta}{2}}
の関係式を使って、
\begin{align} \frac{d\theta}{dt} &=\pm\sqrt{C_0+2\omega^2\left(1-2\sin^2{\frac{\theta}{2}}\right)} \\ &=\pm \sqrt{(C_0+2\omega^2)}\sqrt{1-\frac{4\omega^2}{C_0+2\omega^2}\sin^2\frac{\theta}{2}} \\ &=\pm 2\omega k\sqrt{1-\frac{1}{k^2}\sin^2\frac{\theta}{2}} \end{align}
となります。ここで\displaystyle \frac{1}{k}=\sqrt{\frac{4\omega^2}{C_0+2\omega^2}}という量kを定義しました。
なぜ\frac{1}{k}なの?という理由は後々都合がよくなるからです。深い意味はありません。

運動方程式は
\begin{align} \frac{d\theta}{dt}=\pm 2\omega k\sqrt{1-\frac{1}{k^2}\sin^2\frac{\theta}{2}} \end{align}
とここまで簡略化されました。積分します。
\displaystyle \pm\frac{1}{2\omega k}\int_{\theta_0}^{\theta} \frac{1}{\sqrt{1-\frac{1}{k^2}\sin^2\frac{\theta}{2}}}d\theta=\int_{t_0}^{t} dt
となります。積分区間は初期時刻t_0に角度\theta_0にいた質点が時刻tには角度\thetaにいることを意味します。
通常、積分区間は積分定数C_1を使って表します。この積分区間は以下の変形が可能です。丁寧にやると、
\begin{align} \int_{\theta_0}^{\theta} f(\theta)d\theta &= \int_{t_0}^{t} g(t) dt \\ \int_{\theta_0}^{0} f(\theta)d\theta+\int_{0}^{\theta} f(\theta)d\theta &= \int_{t_0}^{0} g(t) dt+\int_{0}^{t} g(t) dt \\ \int_{0}^{\theta} f(\theta)d\theta &= \int_{0}^{t} g(t) dt + C_1 \end{align}
です。よって
\displaystyle \pm\frac{1}{2\omega k}\int_{0}^{\theta} \frac{1}{\sqrt{1-\frac{1}{k^2}\sin^2\frac{\theta}{2}}}d\theta=\int_{0}^{t} dt + C_1
を解けば良いことになります。左辺の積分を解きましょう。
変数変換\frac{1}{k}\sin\frac{\theta}{2}=\sin{\phi}を行います。

\begin{align} \frac{d}{d\phi}\left(\frac{1}{k}\sin\frac{\theta}{2}\right)&=\frac{d}{d\phi}\sin\phi \\ \frac{1}{k}\cos\left(\frac{\theta}{2}\right)\frac{1}{2}\frac{d \theta}{d\phi} &= \cos\phi \\ \frac{d \theta}{d\phi} &= 2k \frac{\cos{\phi}}{\cos\frac{\theta}{2}} \end{align}
となります。積分区間は
\begin{align} &\theta\ \ | 0 \rightarrow \ \theta \\ &\phi\ \ | 0 \rightarrow \ \phi(=\sin^{-1}(\frac{1}{k}\sin{\frac{\theta}{2}})) \\ \end{align}
と変化します。左辺に代入すれば、
\begin{align} \pm\frac{1}{2\omega k}\int_{0}^{\theta} \frac{1}{\sqrt{1-\frac{1}{k^2}\sin^2\frac{\theta}{2}}}d\theta &=\pm\frac{1}{2\omega k}\int_{0}^{\phi} \frac{1}{\sqrt{1-\sin^2\phi}}2k\frac{\cos{\phi}}{\cos\frac{\theta}{2}} d\phi \\ &=\pm \frac{1}{\omega} \int_{0}^{\phi} \frac{1}{\cos\frac{\theta}{2}}d\phi \\ &=\pm \frac{1}{\omega} \int_{0}^{\phi} \frac{1}{\sqrt{1-k^2\sin^2\phi}}d\phi \\ \end{align}
となります。故に、
\displaystyle \int_{0}^{\phi} \frac{1}{\sqrt{1-k^2\sin^2\phi}}d\phi =\pm\omega t+C_1    …(A)
と表されます。左辺の積分は楕円積分と呼ばれています。

第一種Legendre楕円積分は
\displaystyle F(\phi, k)=\int_0^{\phi} \frac{1}{\sqrt{1-k^2\sin^2\phi}}d\phi
と定義されます。この逆関数のsinはJacobiの楕円関数sn(u,k)をして定義されています。
どういうことかといえば、u=F(\phi, k)のとき、\phi=F^{-1}(u,k)であるので、
\sin\phi=\sin(F^{-1}(u,k))=sn(u,k)
という関数を定義するのです。

すると、今、\pm \omega t+C_1=F(\phi, k)であるため、
\begin{align} \pm \omega t+C_1&=F(\phi, k)=u \\ \phi&=F^{-1}(\pm \omega t+C_1, k) \\ \sin\phi &=\sin(F^{-1}(u,k)) \\ \sin\phi &=sn(u,k) \end{align}
となります。

\sin\phi=\frac{1}{k}\sin\frac{\theta}{2}なので求めたかった\theta(t)
\begin{align} \sin{\frac{\theta}{2}}=k\ sn(\pm\omega t+C_1, k) \\ \theta(t)=2\sin^{-1}(k\ sn(\pm\omega t+C_1, k)) \end{align}
と求められます。

[adsense2]

Jacobiの楕円関数は数値計算的にはランデンの方法という方法で計算されるそうです。
さて、式(A)を見ていると、kが1以上の場合、なんかおかしいことに気が付くかと思います。もしもkが1以上の場合、被積分関数のルートの中が負になることがあり得るんじゃ・・・?
実はこの状況が生まれることはありません。これは積分範囲が今0\sim\phiになっていますが、\phi=\sin^{-1}(\frac{1}{k}\sin{\frac{\theta}{2}})であらわされる区間なので、kの値に応じてとりうる積分範囲は必ずルートの中が負にならない範囲に収まるのです。

ただ、上の問題を心配する必要が無いようにkの範囲は0\le k\le 1に通常は制限されます。
k=2の楕円関数は存在するのに範囲が制限されて、計算どうするの?という疑問が生まれます。
これを解決するのはJacobiの楕円関数snの性質を利用します。
Jacobiの楕円関数は三角関数の拡張であるとみなせるので、その性質は非常に多岐にわたります。
その中の一つにこういう関係式があります。
\displaystyle sn(u,\frac{1}{k})=k\ sn(\frac{u}{k},k)
これを使うとkが1より大きい時でも変換することができます。

具体的には、
\begin{eqnarray} \theta(t)= \left\{ \begin{aligned} & 2\sin^{-1}(k\ sn(\pm\omega t+C_1, k)) \dots (0 \lt k \lt 1) \\ & 2\sin^{-1}(sn(k \omega t+C_1),\frac{1}{k}) \dots (1\lt k) \end{aligned} \right. \end{eqnarray}
となります。運動の状態と対応させるのであれば、一つ目の式(k\lt 1)の場合は往復運動をし、(k\gt 1)の時はぐるんぐるん回る運動に対応します。

1往復、または1回転に要する時間Tは
\begin{eqnarray} T= \left\{ \begin{aligned} &\frac{4}{\omega}K(k) \dots (0 \lt k \lt 1) \\ &\frac{2}{k\omega}K(\frac{1}{k}) \dots (1 \lt k1) \end{aligned} \right. \end{eqnarray}
として与えられます。ここでK(k)は第1種完全楕円積分で、
\displaystyle K(k)=\int_0^1 \frac{1}{\sqrt{(1-x^2)(1-k^2 x^2)}}dx
という定積分です。


例えば、微分方程式
\displaystyle \frac{d^2\theta}{dt^2}=-\sin\theta

t=0,\; \theta=0; \frac{d\theta}{dt}=1.9 (k=0.95)
の初期条件の下解いた場合、
1周期T
T=4K(0.95)=10.360044923498004876778\cdots
となります。
この値はwolfram alphaから求めました。
4*EllipticK[0.95*0.95] wolfram alpha
wolfram alphaでの第1種完全楕円積分は
\displaystyle EllipticK[k’]=\int_0^1 \frac{1}{\sqrt{(1-x^2)(1-k’ x^2)}}dx
として定義されているので注意が必要です。


参考文献


森口繁一, 宇田川銈久, 一松信著『数学公式Ⅲ 特殊関数』岩波書店(1960)p.40~45
Hamilton の正準運動方程式 -epii’s physics notes
KIT -数学ナビゲーション
振り子と楕円関数
Jacobi Elliptic Functions — from Wolfram MathWorld

fortranコンパイラのコマンド

ここではgfortranとifortでのコンパイル時に指定できるコマンド(最適化とデバッグ、mklの使用とopenmp)について記述します。

ここで対象とするファイル名をmain.f90にします。

最適化


  • gfortran
    gfortran -O3 -fopenmp main.f90
  • ifort
    ifort -mkl -openmp -xHOST -ipo -O3 main.f90

としてコンパイルするといいかと。並列計算がいらない場合は-fopenmp, -openmpのオプションを消してください。

デバッグ


デバッグとはバグを取り除いて望み通りのものにする作業です。

  • gfortran
    gfortran -Wall -fbounds-check -O -Wuninitialized -ffpe-trap=invalid,zero,overflow -fbacktrace -g main.f90
  • ifort
    ifort -mkl -check all -warn all -gen_interfaces -fpe0 -ftrapuv -traceback -g main.f90

とするのがいいでしょう。

※gfortranでは、オプション -ffpe-trap=invalid,zero,overflowが悪さすることがあるようです。fortranでのエラーメモにも書いたように、エラーを起こすはずがないコードでもなぜかエラーと言われることがあります。その時は

gfortran -Wall -fbounds-check -O -Wuninitialized -fbacktrace -g main.f90

でデバッグをやると良いでしょう。

[adsense1]

gfortranでMKLを使う


例えばgfortranコンパイラでIntel® Math Kernel Library (Intel® MKL) を並列計算も含めて使いたいとします。

今、こういう状況にあるとします。

MKLがあり、/opt/intel/の中にあるディレクトリ”mkl”
が存在する状況下で
gfortranコンパイラ
がある。

とします。MKLのバージョンによってコンパイルオプションが異なるので注意してください。
Intel® Math Kernel Library Link Line Advisorを見ればどういうコマンドを打てばいいのかを教えてくれます。

ここではMKLのバージョンが
Ver. 11.1

Ver. 10.2.5.035
の場合の具体例を載せます。

MKLver11.1の場合
gfortran -fdefault-integer-8 -fopenmp -m64 -I/opt/intel/mkl/include -Wl,--start-group /opt/intel/mkl/lib/intel64/libmkl_gf_ilp64.a /opt/intel/mkl/lib/intel64/libmkl_core.a /opt/intel/mkl/lib/intel64/libmkl_gnu_thread.a -Wl,--end-group -ldl -lpthread -lm main.f90
MKLver10.2.5.035の場合
gfortran -fopenmp -m64 -I/opt/intel/mkl/include /opt/intel/mkl/lib/em64t/libmkl_solver_lp64.a -Wl,--start-group /opt/intel/mkl/lib/em64t/libmkl_gf_lp64.a /opt/intel/mkl/lib/em64t/libmkl_core.a /opt/intel/mkl/lib/em64t/libmkl_gnu_thread.a -Wl,--end-group -lpthread -lm main.f90

※もしかしたら/opt/intel/mkl/ではなくて、/opt/intel/mkl/10.2.5.035/じゃないとだめかもしれません。

バージョンが違ったらどうするかはIntel® Math Kernel Library Link Line Advisorを参考にしてください。上のものを得るには下のようにすれば大丈夫です。

 Select Intel® product: (MKLのバージョンは?)  Intel(R) MKL10.2
 Select OS: (OSは何?)  Linux*
 Select usage model of Intel® Xeon Phi™ Coprocessor: 
(走らせようとしているコンピュータのcpuの種類がXeon Phi っていうcpuが数百個あるような特別なもの?)
 –
 Select compiler: (何のコンパイラでMKLを使いたいの?, gfortranだったらGNU fortran.)  GNU Fortran
 Select architecture:(良く分かりません!Help me →I don’t use any flag と進んで、”getconf LONG_BIT” と端末で打って何て表示されたかで入力しました。多分cpuが何ビットか、かな?)  Intel(R) 64
 Select dynamic or static linking:

(動的リンク(dynamic)がいい?静的リンク(static)がいい?
    動的→汎用性あり、容量は軽い、動作は遅い
    静的→汎用性なし、容量は重い、動作は早い )
 Static
 Select interface layer:(プログラムの中に整数型で2^31越えてるものはある?ない?越えてなかったら32bit integer, 越えてたら64bit integer)  LP64 (32-bit integer)
 Select sequential or multi-threaded layer:
(mklのサブルーチンを1つのCPUだけ(sequential)で走らせたい?それとも並列(multi-threaded)にする?)
 multi-threaded
 Select OpenMP library:(openMPの種類は何使う?intel か GNUか→僕は両方試したら、intelでは動かず、GNUでは動きました。)  GNU (libgomp)
 Select cluster library: ここから下は何を言っているのか分かりませんでした。空白でもいけたから多分直接は関係ないオプションなのだろうと。  

次に下の方で作られたオプションをgfortranコンパイル時にくっつける。
$(MKLROOT)/libem64t/libmkl_solver_lp64.a -Wl,–start-group $(MKLROOT)/libem64t/libmkl_gf_lp64.a $(MKLROOT)/libem64t/libmkl_core.a $(MKLROOT)/libem64t/libmkl_gnu_thread.a -Wl,–end-group -lpthread -lm

-fopenmp -m64 -I$(MKLROOT)/include
が作られるかと思います。この2つをくっつければokです。
$(MKLROOT)はディレクトリ”mkl”へのパスです。シェルスクリプトの場合、()は要りません。

ここで、この通りコンパイルしようとすると/libem64t/というディレクトリはないよ、っていわれて出来ませんでした。
実際にたどってみると/lib/em64t/ならあって、これに変えたらコンパイルが成功しました。ケースバイケースかな。

参考


Fortranデバッグ用オプション

[adsense2]