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

(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 !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