Fortranで美しいフラクタル構造を描く 2022年4月24日

目次 1. モチベーション 2. フラクタルとは 3. 代表的なフタクタル構造 4. Fortranコードで実装 5. 可視化(VTK) 6. 今後 Fortranへの書き直し VTKFortranライブラリの使用

モチベーション 学生時代にFortranで数値計算をしていた。 最近「C＆Fortran 演習で学ぶ数値計算」が出てプログラミング勉強に火が付いた。

フラクタルとは 自然科学・工学は数式で現象を理解してきた。 しかし、従来のアプローチでは解明できない現象を複雑系と呼ばれる。 流体現象 生体系の現象 経済の現象 •現象自体が理解されていない ⇒数理モデルにすることができない。 •多くの現象の複合 ⇒多くの微分方程式を解く必要があ るため複雑で解けない。 （解いても理解が難しい） 乱流遷移 DANの配列 株価の変動

フラクタル構造 そうは言っても何かパターンがある。 べき乗 𝑓(𝑥) = 𝑎𝑥−𝑡 ｘ 𝑓 𝑎𝑥−𝑡 𝑥を𝑏倍する→𝑏𝑥 ⇔ 𝑓(𝑏𝑥) = 𝑎𝑏−𝑡𝑥−𝑡 = 𝑐𝑥−𝑡 パターン構造 𝑏ｘ 𝑓 𝑐𝑥−𝑡 x軸を拡大 元と同じ形になった。 元と同じ形

実現象でのパターン構造 乱流現象 地震規模と頻度 山崩れ規模と頻度 単語の頻度 動物の斑点模様 リアス式海岸 葉のフラクタル構造 人の介在無しに組織が自己形成をしている。 自己組織化という。 ミクロな視点では、ある決まったルールに従ってマクロな視点ではパターンを形成している。 決まったルールを解明できれば、現象を数理モデルで再現できる。 大小さまざまな渦構造 コルモゴロフ則 http://www.pricewalk.net/concept.html

代表的なフタクタル構造 Fortran Paraview ソルバ ポスト マンデルブロー集合 ジュリア集合 バーンズリーのシダ

「マンデルブロー集合」のアルゴリズム マンデルブロー集合 𝑧𝑛+1 = 𝑧𝑛 2 + 𝐶 初期値𝑧0 で定義される漸化式において𝑛 → ∞の極限で 𝑧𝑛 が発散しない複素数𝐶 = 𝑎 + 𝑏𝑖全体の集合 【例】 𝑧1 = 𝑧0 2 + 𝐶 = 𝐶 𝑧2 = 𝑧1 2 + 𝐶 = 𝑧0 2 + 𝐶 2 + 𝐶 𝑧3 = 𝑧2 2 + 𝐶 = 𝑧1 2 + 𝐶 2 + 𝐶 = 𝑧0 2 + 𝐶 2 + 𝐶 2 + 𝐶 ・・・ 𝑎 𝑏 複素数平面上 発散=0 収束=1 発散=0 収束=1 漸化式

「ジュリア集合」のアルゴリズム ジュリア集合 𝑧𝑛+1 = 𝑧𝑛 2 + 𝐶 𝐶を固定。漸化式において𝑛 → ∞の極限で 𝑧𝑛 が発散 しない複素数初期値𝑧0 = 𝑎 + 𝑏𝑖全体の集合 【例】 𝑧1 = 𝑧0 2 + 𝐶 = 𝐶 𝑧2 = 𝑧1 2 + 𝐶 = 𝑧0 2 + 𝐶 2 + 𝐶 𝑧3 = 𝑧2 2 + 𝐶 = 𝑧1 2 + 𝐶 2 + 𝐶 = 𝑧0 2 + 𝐶 2 + 𝐶 2 + 𝐶 ・・・ 𝑎 𝑏 複素数平面上 発散=0 収束=1 漸化式

「バーンズリーのシダ」のアルゴリズム バーンズリーのシダ 新たな(𝒙, 𝒚)座標 確率1% 確率85% 確率7% 確率7% 𝑥 𝑦 𝑥, 𝑦 = (0,0) 1 2 3 𝒇𝒊 𝒙𝒊 , 𝒚𝒊 = 𝟏. 𝟎 プロットする

Fortranコードで実装 マンデルブロー集合 ジュリア集合 バーンズリーのシダ https://takun- physics.net/13663 https://note.com/kamakiriphysic s/n/n11859282fa79 https://takun- physics.net/13730/

可視化(VTK) !--------------------------------------------------------------------------------- -- !output of movie file if ( mod(m,1000)==0 ) then write(filename,"(a,i5.5,a)") "data/barsley_",int(m),".vtk" open(10,file=filename) write(10,"('# vtk DataFile Version 3.0')") write(10,"('test')") write(10,"('ASCII ')") write(10,"('DATASET STRUCTURED_GRID')") write(10,"('DIMENSIONS ',3(1x,i3))") Nx, Ny, 1 write(10,"('POINTS ',i9,' float')") Nx*Ny*1 do j=0,Ny-1 do i=0,Nx-1 write(10,"(3(f9.4,1x))") x_(i), y_(j), 0.d0 end do end do write(10,"('POINT_DATA ',i9)") Nx*Ny*1 !date input write(10,"('SCALARS barsley float')") write(10,"('LOOKUP_TABLE default')") do j=0,Ny-1 do i=0,Nx-1 write(10,*) f(i,j) end do end do close(10) end if end do # vtk DataFile Version 3.0 test ASCII DATASET STRUCTURED_GRID DIMENSIONS 512 512 1 POINTS 262144 float -3.0000 -3.0000 0.0000 -2.9883 -3.0000 0.0000 -2.9766 -3.0000 0.0000 -2.9648 -3.0000 0.0000 -2.9531 -3.0000 0.0000 -2.9414 -3.0000 0.0000 -2.9297 -3.0000 0.0000 -2.9180 -3.0000 0.0000 -2.9062 -3.0000 0.0000 POINT_DATA 262144 0.0000000000000000 0.0000000000000000 0.0000000000000000 ③ ② ① ① ③ ② ④ ④ ⑥ ⑤ ⑤ ⑥ • https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf • https://itandcfd.com/vtk_file_format/241/ data/ファイル名_0001.vtk

今後  Fortranへの書き直し integer, parameter :: Nmax=500000, Nx=512, Ny=512 double precision, parameter :: Lx=6.d0, Ly=6.d0, dx=Lx/dble(Nx), dy=Ly/dble(Ny) complex(kind(0d0)), parameter :: uso=(0.d0,1.d0) integer :: i, j, n, m, ix, iy real(kind=8)::random_value double precision:: x, y double precision, allocatable, dimension(:) :: x_, y_ double precision, allocatable, dimension(:,:) :: f character(len=40) filename ・・・ do i=0, Nmax call random_seed !乱数の初期値 call random_number(random_value) ! 座標変換 if (0<= random_value .and. random_value<0.01) then x = 0.d0 y = 0.16d0*y ix = int(Nx/2 + x*Nx/10) iy = int(y * Ny /12) f(ix,iy) = 1.0d0  VTKFortranライブラリの使用 module subtoutine, function https://github.com/szaghi/VTKFortran いまさら Fortran？