PROGRAM EXERCISE47 IMPLICIT NONE ! 2次元波動方程式を解くプログラム INTEGER, PARAMETER :: Nx = 100, Ny = 100 ! x方向とy方向の格子点の数 DOUBLE PRECISION, DIMENSION(0:Nx+1, 0:Ny+1) :: u, uold, unew DOUBLE PRECISION :: x(1:Nx), y(1:Ny) INTEGER :: i,j, ti DOUBLE PRECISION :: lx, ly, dx, dy, c DOUBLE PRECISION :: ax, ay, dt, x0, y0, temp CHARACTER(50) :: filename_header, filename, filename_extension, tindex ! x方向、y方向の空間サイズを100(m)とする lx = 100.0d0; ly = 100.0d0 ! x方向とy方向の格子間隔を決める dx = lx / (Nx+1.0d0); dy = ly / (Ny+1.0d0) c = 10.0d0 ! 伝搬速度(m/s) dt = 1.0d-3 ! 時間ステップ(s) ! von Neumannの安定性解析を行って条件が満たされなければエラーにしておく IF( c*c*dt*dt*(1.0/dx**2 + 1.0d0/dy**2) .GT. 1.0d0) THEN WRITE(0,*) "r = ", c*c*dt*dt*(1.0/dx**2 + 1.0d0/dy**2) STOP END IF ! 座標を配列に設定 DO i = 1, Nx x(i) = i*dx END DO DO j = 1, Ny y(j) = j*dy END DO ! 配列を初期化 u = 0.0d0 uold = 0.0d0 unew = 0.0d0 ! 初期関数のパラメターを設定 ax = 3.0d0; ay = 3.0d0 ! 広がりパラメター x0 = 50.0d0; y0 = 50.0d0 ! 中心の座標 ! 初期関数(ガウス型)を設定 DO i = 1, Nx DO j = 1, Ny u(i,j) = EXP( - (x(i)-x0)**2/ax**2 - (y(j)-y0)**2/ay**2) END DO END DO ! 境界条件(自由端) u(0, 0:Ny+1) = u(1, 0:Ny+1) u(Nx+1,0:Ny+1) = u(Nx,0:Ny+1) u(0:Nx+1,0) = u(0:Nx+1, 1) u(0:Nx+1,Ny+1) = u(0:Nx+1,Ny) ! 時刻のカウンタをゼロに設定 ti = 0 ! 初期時刻での関数を出力するファイル名を作成 filename_header = "wave_" ! ファイル名冒頭部分 filename_extension = ".dat" ! ファイル名の拡張子部分 WRITE(tindex,'(I7.7)') ti ! 時刻カウンタを7桁表示して文字型配列tindexに代入 filename = TRIM(filename_header) // TRIM(tindex) // TRIM(filename_extension) !文字列連結してファイル名を作る ! ファイルを開いて初期関数を出力 OPEN(10, FILE=filename, ACTION="WRITE") DO j = 1, Ny DO i = 1, Nx WRITE(10, '(3F18.10)') x(i), y(j), u(i,j) END DO END DO CLOSE(10) DO ! 時刻カウンタを1増やす ti = ti + 1 IF( ti .EQ. 1) THEN unew = u ! t=\Delta tではt=0と同じとする ELSE ! ti > 1 DO j = 1, Ny DO i = 1, Nx temp = (u(i+1,j)-2.0d0*u(i,j)+u(i-1,j))/dx**2 + (u(i,j+1)-2.0d0*u(i,j)+u(i,j-1))/dy**2 unew(i,j) = 2.0d0*u(i,j) - uold(i,j) + dt**2 * c * temp END DO END DO ! 境界条件(自由端) unew(0, 0:Ny+1) = unew(1, 0:Ny+1) unew(Nx+1,0:Ny+1) = unew(Nx,0:Ny+1) unew(0:Nx+1,0) = unew(0:Nx+1,1) unew(0:Nx+1,Ny+1) = unew(0:Nx+1,Ny) END IF ! 100回に一回関数をファイルに書き出す IF( MOD(ti, 100) .EQ. 0) THEN !出力ファイル名を作成 WRITE(tindex,'(I7.7)') ti/100 !ti/100を7桁文字配列に変換 filename = TRIM(filename_header) // TRIM(tindex) // TRIM(filename_extension) !ファイル名を連結して作成 ! ファイルを開いて関数を出力 OPEN(10, FILE=filename, ACTION="WRITE") DO j = 1, Ny DO i = 1, Nx WRITE(10, '(3F18.10)') x(i), y(j), u(i,j) END DO WRITE(10,*) "" ! gnuplot with linesで出力するときはxあるいはyの座標値を変えるときに空行を入れておく(線で結ばれるのを防ぐ) END DO CLOSE(10) PRINT *, "ti = ", ti ! 計算途中経過を出力 END IF ! 次のti計算に向けてuoldとuを更新 uold = u u = unew ! 50000ステップで終了 IF(ti .EQ. 50000) EXIT END DO END PROGRAM EXERCISE47