fortranで台形積分(octave,matlabのcumtrapzを作る)

#はじめに
octaveやmatlabでは、台形積分を

Q=cumtrapz(Y)
Q=cumtrapz(X,Y)

と簡単に書くことができます。Octaveでは、積分結果である配列Qのサイズは事前に定義する必要がありません。fortranでも、allocatable配列を使えば配列サイズを事前に定義せずにこれと同等のことを実現することができます。

#実行可能プログラム
0〜10の区間を等間隔に分割した時の分割点をxとし、xに対応するyの値を0〜1の範囲で生成した乱数とします。このとき、yの値をxで台形積分してみます。
これをoctaveで解く場合、次のように書くことができます。

octaveの場合

x=linspace(0,10,1001);
y=rand(1,1001);
ytrap0=cumtrapz(y);
ytrap1=cumtrapz(x,y);

これと同じ計算をfortranでやる場合、台形積分を自分で書くか、適切な ライブラリを用いることになります。自分で書くなら、octave(やmatlab)と同じような使い勝手の関数を作ることができます。これらの関数を定義したとして、fortranでは次のようになります。

fortranの場合

program testoctave
  use octave    !cumtrapzとlinspaceを定義したモジュールを使う
  implicit none !暗黙の型宣言を無効化
  real(kind=real64),dimension(:),allocatable :: x,y !入力
  real(kind=real64),dimension(:),allocatable :: ytrapz0,ytrapz1,ytrapz2 !出力
  !integer :: i !カウンタ
  !入力データの準備
  x=linspace(0d0,10d0,1001)                    !octave関数と同じもの
  allocate(y(size(x))); call random_number(y)  !積分対象データの生成
  !積分の実施(引数の与え方3パターン)
  ytrapz0=cumtrapz(y)         !yのデータから台形積分(xの間隔は1)
  ytrapz1=cumtrapz(y,x=x)     !x,yのデータペアから台形積分
  ytrapz2=cumtrapz(y,dx=1d-2) !xのデータが等間隔の場合
  !---データ書き出し
  !do i=1,size(x)
  !  write(*,*) x(i),y(i),ytrapz0(i),ytrapz1(i),ytrapz2(i)
  !enddo
end program

データ定義部分がどうしても必要ですが、これを除けばoctaveとfortranでほぼ同じ表現となっています。

#台形積分を実行する関数の定義
linspaceとcumtrapzの定義モジュールを次のように書いてみました。

cumtrapzとlinspaceの定義

module octave
  use iso_fortran_env, only : real64
  implicit none
  public
  contains
    !=============================================
    !台形積分
    pure function cumtrapz(y,x,dx) result(yout)
      real(kind=real64),dimension(:),intent(in) :: y
      real(kind=real64),dimension(:),intent(in),optional :: x
      real(kind=real64),intent(in),optional :: dx
      real(kind=real64),dimension(size(y)) :: yout
      integer :: i,n
      n=size(y)
      yout(1)=0.0_real64
      if(present(x))then
        !x座標値が与えられている場合
        do concurrent(i=2:n)
          yout(i)=(y(i-1)+y(i))*0.5_real64*(x(i)-x(i-1))
        enddo
      elseif(present(dx))then
        !x座標値の間隔がdxで一定となる場合
        do concurrent(i=2:n)
          yout(i)=(y(i-1)+y(i))*0.5_real64*dx
        enddo
      else
        !x座標値は単位値1とする
        do concurrent(i=2:n)
          yout(i)=(y(i-1)+y(i))*0.5_real64
        enddo
      endif
      !次のループには、計算の順序に依存がある
      do i=2,n
        yout(i)=yout(i-1)+yout(i)
      enddo
    end function
    !=============================================
    !x0〜x1区間をnd個のデータで等分割
    pure function linspace(x0,x1,nd) result(x) 
       real(kind=real64),intent(in) :: x0
       real(kind=real64),intent(in) :: x1
       integer,optional,intent(in) :: nd
       real(kind=real64),dimension(:),allocatable :: x
       integer :: n,i,ierr
       real(kind=real64) :: dx 
       if(present(nd))then
         n=nd
         if(n<1) n=1 !負の数はダメ
       else
         n=100 !デフォルト値
       endif
       allocate(x(n))
       if(n==1)then
         x(1)=x1
         return
       endiF
       dx=(x1-x0)/dble(n-1)
       do i=1,n
         x(i)=x0+dx*dble(i-1)
       enddo
     end function
end module

上のように書くと、引数に応じて関数の返り値である配列のサイズは自動的に変わります。さらに、関数の返り値を受け取る側の配列がallocatableであれば、配列サイズはコンパイラが自動的に決定してくれます。

#fortranプログラムの実行結果
x,yのデータは次の通り。

これをcumtrapzを使って積分すると次のようになりました。
(cumtrapz0はxのデータ間隔が1となっているので、これを調整するために0.01を掛けてプロット。)

#おわりに
台形積分をfortranで簡単に実行するため、Octaveと類似の使い勝手となる関数を作りました。このような関数を一度作ってしまえば、メインプログラムを非常に簡単に書くことができます。