program diffusion implicit none real,parameter :: x0=-1.0,x1=1.0 real :: dx,dt integer,parameter :: kosuux=128 integer,parameter :: kosuut=200 integer,parameter :: nint=10 real,dimension(0:kosuux) :: x,ondo real,dimension(0:kosuux) :: ondo2d integer :: i,k,n ! dx = (x1-x0)/real(kosuux) ! x 方向の「小区間」の数 dt = 0.25*dx*dx ! 時間刻み (x方向の刻み幅の2乗に比例させる必要がある) do i = 0,kosuux x(i) = x0+dx*real(i) end do do i = 0,kosuux ! 初期条件を設定 if (abs(x(i)) <= 1.0/32.0 ) then ! |x| <= 1/32 の範囲だけ ondo(i) = 1.0 ! ondo = 1.0 としておく else ondo(i) = 0.0 ! その範囲外では ondo = 0.0 にしておく end if end do k = 0 ! お絵描き用サブルーチンを呼び call oekaki(k,x,ondo) ! 初期条件での ondo の分布を書き出させる do n = 1,kosuut ! 時間発展を計算するループの始まり do i = 1,kosuux-1 ondo2d(i) = & ! 現在の時刻での「温度」の ((ondo(i+1)-ondo(i))/dx & ! 2階微分をとる -(ondo(i)-ondo(i-1))/dx )/dx ! (1階微分をさらに微分する) end do do i = 1,kosuux-1 ! 次の時点での「温度」の値を計算し、 ondo(i) = ondo(i)+dt*ondo2d(i) ! 配列 ondo の中身を更新 end do ! (更新の手順はオイラー法に準じる) ondo(0) = 0.0 ! この2点の値は常に変化しないものとする ondo(kosuux) = 0.0 ! (要するに「境界条件」) if (mod(n,nint) == 0) then ! 時間発展ループが nint 回動くたびに k = n/nint ! お絵描き用サブルーチンを呼び call oekaki(k,x,ondo) ! その時刻での ondo の分布を書き出させる end if ! ちなみに mod は「あまり」を求める関数 end do ! 時間発展を計算するループの終わり stop contains subroutine oekaki(k,x,ondo) integer,intent(in) :: k real,dimension(0:kosuux),intent(in) :: ondo,x character(len=2) :: ck2 ! 長さ2の文字型変数 integer :: i ! 出力用ディレクトリがなければ作成 call system('if [ ! -d diffusion ]; then mkdir -p diffusion; fi') write(ck2,'(i2.2)')k ! 整数 k を文字列に変換 (11 -> 11 とか、3 -> 03 とか) open(10,file="diffusion/"//ck2//".txt") ! 書き出すためのファイルを開く do i = 0,kosuux ! ちなみに // は文字列を連結して新しい文字列を作る命令 write(10,'(12f8.4)')x(i),ondo(i) end do close(10) ! 書き出しが完了したらファイルを閉じる end subroutine oekaki end program diffusion