ロジスティック方程式におけるカオス
カオス
カオスは決定論的方程式の初期値敏感性をさします。ここでは、有名なロジスティック方程式 $x_{n+1} = a x_n (1.0- x_n)$ における初期値依存性を見てみることにします。（ただし$0.0\le x \le 1.0,$ $0\le a \le 4.0$）
Fortran の利用
Fortran には TINY という関数があって、これは最小の正規化数を与えます。また、ある浮動小数点数の最隣接の浮動小数点数を得る NEAREST という関数もあります。これは非正規化数も返してくれます。
さてここでは、最小の正規化数とその両隣の浮動小数点数を出発値として反復を行い、それぞれを R,G,B の三色に割り当ててロジスティック曲線を描いてその差を見てみることにします。もし、差がなければ RGB の和である白色となり、差が出ればそれぞれの色（またはその合成）が現れることになります。計算は倍精度で行うことにします。
実行結果
数値
まず、三つの浮動小数点数を表示してみます。念のため二進表現も表示して、最小ビットのみの差であることを確認しておきます。それぞれの数は約 0.5e-323 だけ離れています。
0.22250738585072009-307 0.22250738585072014-307 0.22250738585072019-307
0000000000001111111111111111111111111111111111111111111111111111
0000000000010000000000000000000000000000000000000000000000000000
0000000000010000000000000000000000000000000000000000000000000001
続行するには何かキーを押してください . . .
画像
次にロジスティック曲線の画像を示します。BMP ファイルを MS-PAINT で PNG に直してうｐしました。横軸は $2.4\le a \le 4.0$ の範囲で、縦軸は $0.0\le x \le 1.0$ となっています。
（右クリック 新しいタブで画像を開く で大きなサイズの画像が見られます。）
図から明らかなように、はじめは白い曲線だったものが、右側のカオスが始まる領域では、色の付いた点になっています。これは、倍精度の最小正規化数とその最隣接の浮動小数点数というわずかな出発数値の差が、カオス領域では大きく拡大していることをしてしています。
初期値の違いで大きく異なった時間発展を見せるというカオス現象が、わずか倍精度の 1ulp の差でも現れることが確かめられました。
また、カオスでは初期値での差が指数関数的に拡大していくとされていますが、単精度と倍精度で比較によって定性的に確認できます。単精度と倍精度では初期値での差が指数部の差で特徴づけられますが、そのファクター倍程度に反復数を増やすだけで同様のカオス的振る舞いが見られます。
プログラム
Intel Fortran v.15 Beta 固有のバグ等に苦しめられて、ちょっと乱れました。
parameterized derived type は、型宣言をした Module 内で、実使用サイズと同じサイズの dummy 変数を確保しておかないと、コンパイラ内部エラーが出ます。またこれを回避した後でも debug mode ではランタイムエラーがでて動かないので release mode で実行する必要があります。
さらに浮動小数点数の演算はデフォルトでは速度重視になっていて非正規化数の演算がうまく行われないので、これを strict モードに指定する必要があります。
module m_bmp
implicit none
type :: t_bmp_file_header
sequence
integer(2) :: bfType = transfer('BM', 0_2, 1) ! BitMap
integer(4) :: bfSize ! file size in bytes
integer(2) :: bfReserved1 = 0 ! always 0
integer(2) :: bfReserved2 = 0 ! always 0
integer(4) :: bfOffBits
end type t_bmp_file_header
!
type :: t_bmp_info_header
sequence
integer(4) :: biSize = Z'28' ! size of bmp_info_header ; 40bytes
integer(4) :: biWidth
integer(4) :: biHeight
integer(2) :: biPlanes = 1 ! always 1
integer(2) :: biBitCount
integer(4) :: biCompression = 0 !0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
integer(4) :: biSizeImage
integer(4) :: biXPelsPerMeter = 3780 ! 96dpi
integer(4) :: biYPelsPerMeter = 3780 ! 96dpi
integer(4) :: biClrUsed = 0
integer(4) :: biClrImportant = 0
end type t_bmp_info_header
!
type :: t_rgb
sequence
character :: b, g, r ! order is b g r
end type t_rgb
!
type :: t_bmp(n_width, n_height)
integer, len:: n_width, n_height
type(t_rgb) :: rgb(n_width, n_height)
contains
procedure :: wr => wr_bmp
procedure :: point => point_bmp
end type
!
! bug of intel fortran v.15beta ? not working in debug mode
! parametrized derived type
type(t_bmp(1280, 720)) :: dummy
!
contains
subroutine wr_bmp(bmp, fn)
class(t_bmp(*, *)), intent(in) :: bmp
character(len = *), intent(in) :: fn
type(t_bmp_file_header) :: bmp_file_header
type(t_bmp_info_header) :: bmp_info_header
associate(nx => size(bmp%rgb, 1), ny => size(bmp%rgb, 2))
bmp_file_header%bfSize = 14 + 40 + 0 + nx * ny * 3
bmp_file_header%bfOffBits = 14 + 40
bmp_info_header%biWidth = nx
bmp_info_header%biHeight = ny
bmp_info_header%biBitCount = 24
bmp_info_header%biSizeImage = nx * ny * 3
end associate
open(9, file = fn//'.bmp', form = 'binary', status = 'unknown')
write(9) bmp_file_header
write(9) bmp_info_header
write(9) bmp%rgb
close(9)
return
end subroutine wr_bmp
!
subroutine point_bmp(bmp, ix, iy, rgb)
class(t_bmp(*, *)), intent(in out) :: bmp
integer, intent(in) :: ix, iy
type(t_rgb), intent(in) :: rgb
integer :: ir, ig, ib
ir = iachar(bmp%rgb(ix, iy)%r) + iachar(rgb%r)
ig = iachar(bmp%rgb(ix, iy)%g) + iachar(rgb%g)
ib = iachar(bmp%rgb(ix, iy)%b) + iachar(rgb%b)
bmp%rgb(ix, iy) = t_rgb(achar(ib), achar(ig), achar(ir))
end subroutine point_bmp
end module m_bmp
program Chaos
use m_bmp
implicit none
integer, parameter :: kd = kind(1.0d0)
type (t_rgb), parameter :: blue = t_rgb(achar(255), achar(0), achar(0))
type (t_rgb), parameter :: green = t_rgb(achar(0), achar(255), achar(0))
type (t_rgb), parameter :: red = t_rgb(achar(0), achar(0), achar(255))
real(kd) :: a, a0, a1, x0, x(3)
integer :: i, k, ix(3)
type (t_bmp(1280, 720)) :: fig
!
fig%rgb = t_rgb(achar(0),achar(0),achar(0))
x0 = tiny(0.0_kd)
a0 = 2.4_kd
a1 = 4.0_kd
do k = 0, fig%n_width - 1
a = (a1 - a0) / real(fig%n_width - 1, kd) * k + a0
x = [nearest(x0, -1.0_kd), x0, nearest(x0, +1.0_kd)]
if (k == 0) print '(3e25.17/3(b66.64/))', x, x
do i = 1, 1000
x = a * x * (1.0_kd - x)
end do
do i = 1, 100
x = a * x * (1.0_kd - x)
ix = x * (fig%n_height - 1) + 1
call fig%point(k, ix(1), green)
call fig%point(k, ix(2), red )
call fig%point(k, ix(3), blue )
end do
end do
call fig%wr('chaos')
end program Chaos
補足： rbg のところで、achar や iachar を使っている理由。
Fortran には符号なし整数がないために、1byte 整数型を用いると 128 以上の数が負数になってしまい色々面倒なことが起きます。これを避けるために、文字コードを符号なし 1byte 整数の代用として利用しています。