Julia 集合

複素平面上の点を漸化式の初期値として、漸化式の極値で初期値を分類したものをJulia集合といいます。

ここでは三次方程式 $f(x)=x^3-1=0$ をニュートン法を用いて反復的に解くことを考えます。ニュートン法による反復は漸化式であらわせば、$x_{n+1}=x_n-f(x_n)/f'(x_n)$ すなわち $x_{n+1}=x_n-{{x_n^3-1}\over 3x_n^2}$ となります。

一方、この方程式の解は、$x=1,{{-1\pm\sqrt{3}}\over2}$ の三つと自明的に分かります。

さて初期値として複素平面上の点を与えると、その初期値に応じてニュートン法で得られる解が定まります。三つの解のうちどの値に収束してゆくかを初期値の座標にプロットしてゆくと、フラクタル的な自己相似図形が得られます。

いま三つの解は等価であるので、複素平面上に三回対称な図形が得られます。キャラクターでプロットすると、以下のようになります。

Julia 集合　キャラクター版

ここで、1は$1$、2は${{-1-\sqrt{3}}\over2}$、3は${{-1+\sqrt{3}}\over2}$に収束したことを意味します。4は与えた反復数の範囲ではまだどの解に収束してゆくのか判断できないことを意味しています。漸化式の反復数は最大15回としています。

3333333333333333333333333333333333333333333333333333341111142441111111111111111
3333333333333333333333333333333333333333333333333333343332222333411111111111111
3333333333333333333333333333333333333333333333333333412222223333111111111111111
3333333333333333333333333333333333333333333333333333422222222422111111111111111
3333333333333333333333333333333333333333333333333333322222243411111111111111111
3333333333333333333333333333333333333333333333333333322231111111111111111111111
3333333333333333333333333333333333333333333333321111234111111111111111111111111
3333333333333333333333333333333333333333333333312222233311111111111111111111111
3333333333333333333333333333333333333333333333122222231411111111111111111111111
3333333333333333333333333333333333333333333333122222331111111111111111111111111
3333333333333333333333333333333333333333333333322111111111111111111111111111111
3333333333333333333333333333333333333333321112231111111111111111111111111111111
3333333333333333333333333333333333333331122222433211111111111111111111111111111
3333333333333333333333333333333333333331222222242111111111111111111111111111111
1333333333333333333333333333333333333332222222111111111111111111111111111111111
3114123333342224223333332222421221333332222231111111111111111111111111111111111
1111111113322111111112342311111111111332211111111111111111111111111111111111111
1111111112233111111113243211111111111223311111111111111111111111111111111111111
2114132222243334332222223333431331222223333321111111111111111111111111111111111
1222222222222222222222222222222222222223333333111111111111111111111111111111111
2222222222222222222222222222222222222221333333343111111111111111111111111111111
2222222222222222222222222222222222222221133333422311111111111111111111111111111
2222222222222222222222222222222222222222231113321111111111111111111111111111111
2222222222222222222222222222222222222222222222233111111111111111111111111111111
2222222222222222222222222222222222222222222222133333221111111111111111111111111
2222222222222222222222222222222222222222222222133333321411111111111111111111111
2222222222222222222222222222222222222222222222213333322211111111111111111111111
2222222222222222222222222222222222222222222222231111324111111111111111111111111
2222222222222222222222222222222222222222222222222222233321111111111111111111111
2222222222222222222222222222222222222222222222222222233333342411111111111111111
2222222222222222222222222222222222222222222222222222433333333433111111111111111
2222222222222222222222222222222222222222222222222222413333332222111111111111111
2222222222222222222222222222222222222222222222222222242223333222411111111111111
2222222222222222222222222222222222222222222222222222241111143441111111111111111
2222222222222222222222222222222222222222222222222222222134344444111111111111111
続行するには何かキーを押してください . . .

プログラム

Fortran2003 で配列の index を返す関数は、minloc, maxloc 位しかないのですが、F2008 では findloc という関数が加わります。ここでは適当な専用の関数を定義しました。

    program julia
      implicit none
      integer, parameter :: nx = 79, ny = 35  ! nx should be a multiple of 4 ! 24bit bmp
      complex :: c(nx, ny)

      block    
        integer :: ix, iy
        forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(4.0 / nx * ix - 2.0, 4.0 / ny * iy - 2.0)
      end block

      print '(79i1)', ijulia(c)

      stop      
    contains
      pure elemental integer function ijulia(c) 
        complex, intent(in) :: c
        complex, parameter :: v(3) = [(1.0, 0.0), cmplx(-1.0, sqrt(3.0)) / 2.0, cmplx(-1.0, -sqrt(3.0)) / 2.0]
        complex :: z
        integer :: i, k
        z = c
        do i = 1, 15
          ijulia = locq( abs(z - v) < 0.05 )
          if (ijulia <= 3) exit
          z = z - (z**3 - 1.0) / (3 * z**2)
        end do    
      end function ijulia

      pure integer function locq(q)
        logical, intent(in) :: q(:)
        do locq = 1, size(q)
          if (q(locq)) exit
        end do
      end function locq
    end program julia  

Julia 集合 BMP 版

プログラム

    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
        type(t_rgb), allocatable :: rgb(:, :)
      contains 
        procedure :: wr => wr_bmp
      end type
    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       ! nx shouold be multiles of 4
          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
 ! convert to t_RGB    
      pure elemental type(t_rgb) function to_rgb(ir, ig, ib)
        integer, intent(in) :: ir, ig, ib
        to_rgb = t_rgb(achar(ib), achar(ig), achar(ir))
      end function to_rgb  
    end module m_bmp

    program julia
      use m_bmp
      implicit none
      integer, parameter :: nx = 1280, ny = 1280  ! nx should be a multiple of 4 ! 24bit bmp
      complex :: c(nx, ny)

      block    
        integer :: ix, iy
        forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(4.0 / nx * ix - 2.0, 4.0 / nx * iy - 2.0)
      end block

      block
        type(t_bmp) :: bmp
        bmp%rgb = color( ijulia(c) )
        call bmp%wr('Julia')
      end block    

    contains
      pure elemental integer function ijulia(c) 
        complex, intent(in) :: c
        complex, parameter :: v(3) = [(1.0, 0.0), cmplx(-1.0, sqrt(3.0)) / 2.0, cmplx(-1.0, -sqrt(3.0)) / 2.0]
        complex :: z
        integer :: i, k
        z = c
        do i = 1, 15
          ijulia = locq( abs(z - v) < 0.05 )
          if (ijulia <= 3) exit
          z = z - (z**3 - 1.0) / (3 * z**2)
        end do    
      end function ijulia

      pure integer function locq(q)
        logical, intent(in) :: q(:)
        do locq = 1, size(q)
          if (q(locq)) exit
        end do
      end function locq

      pure elemental type(t_rgb) function color(i)
        integer, intent(in) :: i
        character, parameter :: aff = achar(255), a00 = achar(0)
        type(t_rgb), parameter :: col(4) = [t_rgb(aff, a00, a00), t_rgb(a00, aff, a00), t_rgb(a00, a00, aff), t_rgb(a00, a00, a00)]
        color = col(i)
      end function color

    end program julia  

参考文献：
奥村晴彦 著「C言語による最新アルゴリズム事典 」 技術評論社　(1991)