カオスの初期値依存性

#ロジスティック方程式におけるカオス
###カオス
カオスは決定論的方程式の初期値敏感性をさします。ここでは、有名なロジスティック方程式 $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$ となっています。
![chaos.png](https://qiita-image-store.s3.amazonaws.com/0/44522/82d3fa2e-3bf5-b892-694d-5082559628b7.png)
（右クリック　新しいタブで画像を開く　で大きなサイズの画像が見られます。）

図から明らかなように、はじめは白い曲線だったものが、右側のカオスが始まる領域では、色の付いた点になっています。これは、倍精度の最小正規化数とその最隣接の浮動小数点数というわずかな出発数値の差が、カオス領域では大きく拡大していることをしてしています。

初期値の違いで大きく異なった時間発展を見せるというカオス現象が、わずか倍精度の 1ulp の差でも現れることが確かめられました。

また、カオスでは初期値での差が指数関数的に拡大していくとされていますが、単精度と倍精度で比較によって定性的に確認できます。単精度と倍精度では初期値での差が指数部の差で特徴づけられますが、そのファクター倍程度に反復数を増やすだけで同様のカオス的振る舞いが見られます。

###プログラム
Intel Fortran v.15 Beta 固有のバグ等に苦しめられて、ちょっと乱れました。

parameterized derived type は、型宣言をした Module 内で、実使用サイズと同じサイズの dummy 変数を確保しておかないと、コンパイラ内部エラーが出ます。またこれを回避した後でも debug mode ではランタイムエラーがでて動かないので release mode で実行する必要があります。

さらに浮動小数点数の演算はデフォルトでは速度重視になっていて非正規化数の演算がうまく行われないので、これを strict モードに指定する必要があります。

```fortran
    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
```

```fortran    
    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 整数の代用として利用しています。