3 プログラム

今回は Fortran と C++ と両方出してみよう。初期条件は $x=0.5$ のところ にだけ値がある $\delta$ 関数的なものということにする。

なお、面倒なので $D=1$ ということのする。時間の単位が $1/D$ と思えば同 じことである。

pt

c
c     program parabolic1
c
c     Copyright 2000 J. Makino

      subroutine initparam(u, nx, dx, dt, niter, gmode)
      real*8 u(*), dx, dt, tmax
      integer nx, niter, gmode, i
      write(*,*)'Enter nx, tmax, dt:'
      read(5,*) nx, tmax, dt
      dx = 1d0/nx
      niter = (tmax+dt/2)/dt
      write(*,*)'Enter graphics mode 0: no G'
      write(*,*)'                    1: animation'
      write(*,*)'                    2: no-erase)'
      read(5,*) gmode
      write(6,*) 'nx = ', nx, ' dt = ', dt, ' niter = ', niter
      write(6,*) 'gmode = ', gmode
      do i = 1, nx+1
         u(i) = 0
      enddo
      u((nx+2)/2) = nx
      end

      subroutine push_system(u, unew, nx, dx, dt)
      real*8 u(*), unew(*), dx, dt, lambda
      integer nx, i
      lambda = dt/(dx*dx)
      do i = 2, nx
         unew(i) = u(i) + lambda*(u(i-1)-2*u(i)+u(i+1))
      enddo
      do i = 2, nx
         u(i) = unew(i)
      enddo
      end

      subroutine initgraph
      integer pgopen
      if(pgopen('?') .LE. 0) stop
      call pgask(.false.)
      end
      
      subroutine show_graphics(u, nx, dx, gmode,
     $     xmin, xmax, ymin,ymax, first)
      real*8 u(*), dx
      real*8 xmin, xmax, ymin, ymax
      integer nx, gmode, i, first
      if (gmode .eq. 0) return
      call pgbbuf
      if (first .eq. 1) then
         call pgpage
         call pgenv(real(xmin), real(xmax), real(ymin), real(ymax),0,-2)
      endif
      if (gmode .eq. 1) then
         call pgeras
      endif
      call pgbox('bcnst', 0.0, 0, 'bcnst', 0.0, 0)
      call pglab('X', 'Y', ' ')
      call pgmove(0.0, real(u(1)))
      do i = 2, nx+1
         call pgdraw(real((i-1)*dx), real(u(i)))
      enddo
      call pgebuf
      end
      
      program parabolic
      integer nmax
      parameter (nmax=10001)
      real*8 u(nmax), unew(nmax), dx, dt
      integer nx, niter, gmode, iter
      call initparam(u, nx, dx, dt, niter, gmode)
      if (gmode .ne. 0) call initgraph
      call show_graphics(u, nx, dx, gmode,
     $        0d0, 1d0, 0d0, 10d0, 1)
      do iter = 1, niter
         call push_system(u, unew, nx, dx, dt)
         call show_graphics(u, nx, dx, gmode,
     $        0d0, 1d0, 0d0, 10d0, 0)
      enddo
      end

これを適当なファイルに落して、コンパイル・実行する。コンパイルには

        f77  -o parabolic1  parabolic1.f -L/usr/local/pgplot \
                -lpgplot -L/usr/X11R6/lib -lX11

とすると、 parabolic1 という名前の実行ファイルができる。

ここでは、 PGPLOT というわりと広く使われているグラフィックライブラリを 使う。これは画面、 PSファイル等様々な出力形態を同じプログラムで切替え られるので割合に便利である。

中身の解説は、見ての通りという気もする。次に C++ のプログラムをつける