'************************************************************************
'*  Solve Laplace Equation by relaxation Method: d2T/dx2 + d2T/dy2 = 0  *
'* -------------------------------------------------------------------- *
'* Example: Determine the temperatures T(x,y) in a square plate [x,y]   *
'*          of side length C with following limit conditions:           *
'*          1) T = 0 for sides x = 0 and y = 0.                         *
'*          2) T = x^3 for side y = C                                   *
'*          3) T = 16*x for side x = C  (see sketch below).             *
'*          Other parameters:                                           *
'*          square length C = 4 (any unit)                              *
'*          Number of subdivisions NM = 8 (9 knots per line from 0 to 8)*
'*          integration steps dx or dy = C/NM                           *
'*          weight coefficient w = 1.45                                 *
'*          Maximum residual error ER = 0.01                            *
'*          Maximum number of iterations IT = 40                        *
'* -------------------------------------------------------------------- *
'* SAMPLE RUN:                                                          *
'*                                                                      *
'* Niter  Max. Residual     W                                           *
'*  18      0.0031         1.45                                         *
'* Temperature                                                          *
'*  0.000  0.125  1.000  3.375   8.000  15.625  27.000  42.875  64.000  *
'*  0.000  1.401  3.401  6.586  11.515  18.693  28.509  41.106  56.000  *
'*  0.000  2.078  4.616  8.052  12.783  19.122  27.237  37.040  48.000  *
'*  0.000  2.296  4.932  8.224  12.442  17.774  24.278  31.818  40.000  *
'*  0.000  2.175  4.592  7.469  10.987  15.254  20.283  25.953  32.000  *
'*  0.000  1.814  3.790  6.075   8.782  11.974  15.646  19.711  24.000  *
'*  0.000  1.291  2.681  4.259   6.092   8.213  10.616  13.244  16.000  *
'*  0.000  0.668  1.384  2.189   3.113   4.171   5.359   6.651   8.000  *
'*  0.000  0.000  0.000  0.000   0.000   0.000   0.000   0.000   0.000  *
'*                                                                      *
'* -------------------------------------------------------------------- *
'* REFERENCE:  "Méthode de calcul numérique- Tome 2 - Programmes en     *
'*              Basic et en Pascal By Claude Nowakowski, Edition du     *
'*              P.S.I., 1984" [BIBLI 04].                               *
'*                                                                      *
'*                                 Basic Release By J-P Moreau, Paris.  *
'*                                          (www.jpmoreau.fr)           *
'************************************************************************
'PROGRAM LAPLACE

OPTION BASE 0
                   '    y   T=x^3
IT = 40            '     ------------           
ER = .01           '    |     C      |
OM = 1.45          '    |            |
NM = 8             ' T=0|           C|T=16x
IC = 4             '    |            |
                   '    |            |
DIM T(NM, NM)      '     ------------ x
                   '         T=0
  CLS

  dx = (1! * IC) / NM
  T = 0!
  FOR i = 0 TO NM
    T(i, NM) = (dx * i) ^ 3
    T(NM, i) = 16! * dx * i
  NEXT i
  n1 = NM - 1: l = 0
  PRINT
  PRINT " Niter      Error"
10  rm = 0!
    l = l + 1

    FOR i = 1 TO n1
      FOR j = 1 TO n1
        GOSUB 1000  'call Schema(r,i,j,T)
        T(i, j) = T(i, j) + .25 * OM * r
        IF ABS(r) > rm THEN rm = ABS(r)
      NEXT j
    NEXT i
  IF rm > ER AND l <= IT THEN GOTO 10
  CLS
  PRINT
  IF l < IT THEN
    PRINT " Niter     Max. Residual     W"
    PRINT USING "####      ####.####      ####.##"; l; rm; OM
    PRINT " Temperature"
    FOR j = NM TO 0 STEP -1
      FOR i = 0 TO NM
        PRINT USING "##.###"; TAB(9 * i); T(i, j);
      NEXT i
      PRINT
    NEXT j
  ELSE
    PRINT " Convergence, RMAXI = "; rm
  END IF

END

1000 'calculate residual value r at current point (i,j)
  r = T(i + 1, j) + T(i - 1, j) + T(i, j + 1) + T(i, j - 1) - 4 * T(i, j)
RETURN

'end of file Laplace.bas