More about the Explicit Method

 

$$T(x_0, y_0, t_0 + \Delta t)$$

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

$$\qquad + T(x_0, y_0 + h, t_0) + T(x_0, y_0 - h, t_0)$$

$$\qquad - 4 T(x_0, y_0, t_0) \big) + T(x_0, y_0, t_0)$$

$$\quad = \frac{\alpha}{4} \big(T(x_0 + h, y_0, t_0) + T(x_0 - h, y_0, t_0)$$

$$\qquad + T(x_0, y_0 + h, t_0) + T(x_0, y_0 - h, t_0)$$

$$\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$:

$$\Delta t = \frac{h^2}{4 D}$$ (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.

• There is a thoroughly brain-damaged method called Successive Over Relaxation method (SOR, for short).  
• This method says: make $\alpha > 1$ - this is like increasing $D\Delta t$ so you'll end up converging faster, because your time step and your diffusion constant effectively increase (assuming that you don't change h, i.e., your resolution).
• SOR violates the Courant-Lewy condition and thus must go unstable, eventually, although you may be able to SOR-time-step for a while without a major disaster.
• It is seldom worth the bother.

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