Next: Spectral Methods Up: Fields and Rasters Previous: Solving the Time Dependent

### More about the Explicit Method

 (3.22)

where . Observe that when we get the Jacobi iteration scheme. And, clearly, we can always choose so that :

 (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 - this is like increasing 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 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
!hpfalign raster(:,:) with field !hpf distribute (*, block) :: field

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

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
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: Spectral Methods Up: Fields and Rasters Previous: Solving the Time Dependent
Zdzislaw Meglicki
2001-02-26