計算機基礎B - 10. 関数とグラフの作図

目次

リンク

簡単な例

次のように定義される関数 esin(x) を Fortran の関数として作成し，関数の概形を Gnuplot で図示してみましょう．

esin(x) = sin(x) * exp(-x*x/50)

これを関数副プログラムにすると次のようになります． 結果を表示する write 文の中の '(2E20.12)' という文字列は， 2つの実数型変数を「小数部分12桁・全体20桁」で表示するための「書式 (format)」です． あとで詳しく解説しますので，ここではおまじないと思って使ってください．

! esin(x) を計算する関数
real(8) function esin(x)
  implicit none
  real(8) :: x

  esin = sin(x) * exp(-x*x/50.0d+0)

  return
end function esin

! [-10, 10] の範囲で esin(x) を計算し，結果を出力するメインプログラム
program main
  implicit none
  real(8) :: esin               ! ユーザ定義関数には型宣言が必要
  real(8) :: x
  integer :: i

  do i = -100, 100
     x = i / 10.0d+0
     write(*,'(2E20.12)') x, esin(x)
  end do

  stop
end program main

ソースプログラム（ファイル名を esin.f90 とします）をコンパイルして実行し，データを esin.out に保存します． 以下はターミナルから入力するコマンドの例です． % はシェルのプロンプトなので入力する必要はありません．

% gfortran esin.f90
% ./a.out > esin.out

Gnuplot を使って esin.out をグラフ表示してみましょう． gnuplot> は Gnuplot のプロンプトです． plot 以下のコマンドを入力しましょう．

% gnuplot
gnuplot> plot 'esin.out' w l t 'esin(x)'
gnuplot> quit

w は with の略，l (エル) は line の略，t は title の略です． quit コマンドで Gnuplot を終了します．

小数部分を返す関数

関数 afrac(x) を実数 x の小数部分を返す関数と定義します． 例えば，

afrac(1.2) = 0.2
afrac(2.8) = 0.8
afrac(-3.4) = -0.4

などです． afrac(x) を Fortran の関数副プログラムとして書くには， x の小数点以下を切り捨てる組み込み関数 aint(x) を使うとよいでしょう． aint(x) は次のような計算を行います．

aint(1.2) = 1.0
aint(9.8) = 9.0

以下のプログラムを完成させ，afrac(x) を [-10, 10] の範囲で図示してみましょう．

! 小数部分を返す関数
real(8) function afrac(x)
  implicit none
  real(8) :: x

  !*****************************************************************
  ! この部分を自分で作ってみましょう
  !*****************************************************************

  return
end function afrac

! データ作成のためのメインプログラム
program main
  implicit none
  real(8) :: afrac
  real(8) :: x
  integer :: i

  do i = -100, 100
     x = i / 10.0d+0
     write(*,'(2E20.12)') x, afrac(x)
  end do

  stop
end program main

プログラムが完成したらコンパイルして実行します． 正しく動作するようなら，出力データをファイルに保存し，Gnuplot で図示します． やり方は上の例を参考にしてください．

テイラー展開を計算する関数

指数関数 exp(x) の原点を中心とする n 次までのテイラー展開は次のようになります．

exp(x) = 1 + x + x**2/2! + x**3/3! + ... + x**n/n!

これを計算する関数副プログラム et(n, x) を作ってみましょう． ここで，x**n は x の n 乗，n! は n の階乗を表します． c(k) = x**k / k! と定義すると，

et(n, x) = c(0) + c(1) + c(2) + ... + c(n)

と書くことができます． c(k) は次の漸化式を満たすので，簡単に計算することができます．

c(0) = 1
c(k) = c(k-1) * (x / k)       (k = 1, 2, 3, ...)

do ループの中で c(0), c(1), ..., c(n) を順次計算し，足し合わせていけば，et(n, x) を計算することができます． この部分を書き足し，次のプログラムを完成させてください．

! exp(n) のテイラー展開を n 次まで計算する関数
real(8) function et(n, x)
  implicit none
  integer :: n
  real(8) :: x
  ! ローカル変数
  integer :: k
  real(8) :: cn

  ! et = c(0) + c(1) + ... + c(n) を計算する
  cn = 1.0d+0                   ! c(0) = 1
  et = cn                       ! et = c(0)
  do k = 1, n
     !*****************************************************************
     ! この部分を自分で作ってみましょう
     !*****************************************************************
  end do

  return
end function et

! [-3, 3] の範囲で関数をプロットする
program main
  implicit none
  real(8) :: et
  real(8) :: x
  integer :: i

  do i = -30, 30
     x = i / 10.0d+0
     write(*,'(5E20.12)') x, exp(x), et(0, x), et(5, x), et(10, x)
  end do

  stop
end program main

プログラムをコンパイルして実行すると，各行に x, exp(x), et(0, x), et(5, x), et(10, x) の値が表示されます． exp(x) と et(5, x) のグラフを比較するには次のようにします．

% gfortran et.f90
% ./a.out > et.out
% gnuplot
gnuplot> plot 'et.out' using 1:2 w l, 'et.out' using 1:4 w l

using 1:2 は1番目のデータを x 座標の値，2番目のデータを y 座標の値としてグラフをプロットするという意味です． using 1:4 では1番目と4番目のデータを使います．

同じようにして，exp(x) と et(10, x) のグラフも比べてみましょう．

ルジャンドル多項式

n 次のルジャンドル多項式 Pn(x) を計算する Fortran の関数 lp(n, x) を作成しグラフを描いてみましょう（関数名はエルピーです）． Pn(x) は [-1, 1] の実数 x に対して次の漸化式で定義されます．

P0(x) = 1
P1(x) = x
Pn(x) = [ (2*n-1) * x * Pn-1(x) - (n-1) * Pn-2(x) ] / n      (n = 2, 3, ...)

以下は関数副プログラム lp(n, x) の例です． 不足している部分を書き足してプログラムを完成させてください． n = 0 または 1 の場合，漸化式の計算は不要なので，lp の値をセットしたらすぐに return を実行し，副プログラムから脱出します． n >= 2 の場合は do ループを使って Pn まで漸化式を計算します．

!------------------------------------------------------------
! n 次のルジャンドル多項式 Pn(x)
!------------------------------------------------------------
real(8) function lp(n, x)
  implicit none
  integer, intent(IN) :: n
  real(8), intent(IN) :: x
  ! local
  integer :: k
  real(8) :: pa, pb             ! 漸化式用の変数

  ! n = 0 の場合
  if (n == 0) then
     lp = 1                     ! lp に値をセット
     return                     ! 関数から抜ける
  end if

  ! n = 1 の場合
  if (n == 1) then
     !*****************************************************************
     ! この部分を自分で作ってみましょう
     !*****************************************************************
  end if

  ! n = 2, 3, ... の場合
  pa = 1.0d+0                                     ! P0
  pb = x                                          ! P1
  do k = 2, n
     !*****************************************************************
     ! この部分を自分で作ってみましょう
     !*****************************************************************
  end do

  return
end function lp

!------------------------------------------------------------
! メインプログラム
!------------------------------------------------------------
program main
  implicit none
  integer :: n
  real(8) :: x, lp              ! 関数 lp() は real(8) 型

  do n = -20, 20
     x = n / 20.0d+0
     write(*,'(5E20.12)') x, lp(1, x), lp(2, x), lp(3, x), lp(4, x)
  end do

  stop
end program main

プログラムを実行すると，各行に x, P1(x), P2(x), P3(x), P4(x) の値が出力されます． 結果を lp.out に保存して Gnuplot で表示してみましょう． ターミナルから入力するコマンドは次のようになります．

% gfortran lp.f90
% ./a.tou > lp.out
% gnuplot
gnuplot> plot 'lp.out' using 1:2 w l, 'lp.out' using 1:3 w l,  (実際は改行しない)
              'lp.out' using 1:4 w l, 'lp.out' using 1:5 w l

次回までの課題

n 次のエルミート多項式 Hn(x) を計算する関数副プログラム hp(n, x) を作成せよ． 任意の実数 x に対して Hn(x) は次の漸化式で定義される．

H0(x) = 1
H1(x) = 2 * x
Hn(x) = 2 * x * Hn-1(x) - 2 * (n-1) * Hn-2(x)     (n = 2, 3, 4, ...)

さらに，Gnuplot を用いて [-10, 10] の範囲で関数 Hn(x) * exp(-x*x/2) (n = 0, 1, 2, 3) をプロットせよ．

完成したら高木または TA 田中さんのチェックを受けること．