next up previous index
Next: Spectral Methods Up: Fields and Rasters Previous: Solving the Time Dependent

More about the Explicit Method


 
    $\displaystyle T(x_0, y_0, t_0 + \Delta t)$  
    $\displaystyle \quad =
\frac{D \Delta t}{h^2}
\big(T(x_0 + h, y_0, t_0) + T(x_0 - h, y_0, t_0)$  
    $\displaystyle \qquad + T(x_0, y_0 + h, t_0) + T(x_0, y_0 - h, t_0)$  
    $\displaystyle \qquad - 4 T(x_0, y_0, t_0) \big) + T(x_0, y_0, t_0)$  
    $\displaystyle \quad =
\frac{\alpha}{4}
\big(T(x_0 + h, y_0, t_0) + T(x_0 - h, y_0, t_0)$  
    $\displaystyle \qquad + T(x_0, y_0 + h, t_0) + T(x_0, y_0 - h, t_0)$  
    $\displaystyle \qquad - 4 T(x_0, y_0, t_0) \big) + T(x_0, y_0, t_0)$ (3.22)

where $\alpha = 4 D \Delta t/h^2$. Observe that when $\alpha = 1$ we get the Jacobi iteration scheme. And, clearly, we can always choose $\Delta t$ so that $\alpha = 1$:

\begin{displaymath}\Delta t = \frac{h^2}{4 D}
\end{displaymath} (3.23)

This is, in fact, Courant-Lewy   condition for parabolic equations. You can time-step slower, but you can't time-step faster without going unstable.

The following SOR program lets you experiment with various values of $\alpha$ and various numbers of time steps. Images are dumped every second time step so that you can observe the onset of instability. Observe the use of safety breaks, i.e., the variable field is forced to stay within a certain range between floor and ceiling.

program jacobi

  use hdf
  use dffunc

  implicit none

!hpf$ nosequence

  integer, parameter :: n = 200
  real, parameter :: floor = 30.0, ceiling = 250.0
  real, dimension(n), parameter :: north_boundary = 30.0, &
       east_boundary = 160.0, west_boundary = 160.0, south_boundary = 250.0
  real, dimension(n, n) :: field = 30.0
  character(len=1), dimension(n, n) :: raster
  logical, dimension(n, n) :: mask = .true.
  real :: alpha
  integer :: i, j, status, iterations

!hpf$ align mask(:,:) with field
!hpf$ align raster(:,:) with field
!hpf$ distribute (*, block) :: field

  field(ubound(field, dim=1), :) = east_boundary
  mask(ubound(mask, dim=1), :) = .false.
  field(lbound(field, dim=1), :) = west_boundary
  mask(lbound(mask, dim=1), :) = .false.
  field(:, ubound(field, dim=2)) = north_boundary
  mask(:, ubound(mask, dim=2)) = .false.
  field(:, lbound(field, dim=2)) = south_boundary
  mask(:, lbound(mask, dim=2)) = .false.

  raster=char(int(field))
  status=d8pimg('sor-movie.hdf', raster, n, n, COMP_RLE)

  write(*, '(1a, $)') 'alpha > '; read(*, *) alpha
  write(*, *) 'alpha = ', alpha
  write(*, '(1a, $)') 'iter  > '; read(*, *) iterations
  write(*, *) 'iter  = ', iterations

  do i = 1, iterations
     where (mask)
        field = field + &
             (eoshift(field, 1, dim=1) + eoshift(field, -1, dim=1) &
             + eoshift(field, 1, dim=2) + eoshift(field, -1, dim=2) &
             - 4 * field) * 0.25 * alpha
     end where
     where (field .gt. ceiling) field = ceiling
     where (field .lt. floor) field = floor
     if (mod(i, 2) .eq. 0) then
        raster = char(int(field))
        status = d8aimg('sor-movie.hdf', raster, n, n, comp_rle)
        write(*,*) i, ' -> ', (ichar(raster(j, 10)), j = 1, 5), &
             (ichar(raster(j, 10)), j=96, 105), &
             (ichar(raster(j, 10)), j=196, 200)
     end if
  end do

end program jacobi


next up previous index
Next: Spectral Methods Up: Fields and Rasters Previous: Solving the Time Dependent
Zdzislaw Meglicki
2001-02-26