'**************************************************************************
'* Solve Laplace Equation by relaxation Method: d2T/dx2 + d2T/dy2 = 0 (3) *
'* ---------------------------------------------------------------------- *
'* Example #3: Determine the temperatures T(x,y) in a domain [x,y] made   *
'*             of one square of length C with the following conditions:   *
'*             along OC: T=0, along AB: T=1, along OA and CB: dT/dy = 0   *
'*             (see sketch below).                                        *
'*             Main parameters:                                           *
'*             C: length of square domain.                                *
'*             Number of subdivisions 8 in ox and 8 in Oy.                *
'*             integration steps dx=C/8.                                  *
'*             weight coefficient w = 1.3                                 *
'*             Maximum residual error ER = 0.001                          *
'*             Maximum number of iterations IT = 100                      *
'* ---------------------------------------------------------------------- *
'* SAMPLE RUN:                                                            *
'*                                                                        *
'* Niter  Max. Residual     W                                             *
'*  43      0.0009         1.30                                           *
'* Temperature                                                            *
'*  0.000  0.125  0.249  0.374  0.499  0.624  0.749  0.875  1.000         *
'*  0.000  0.125  0.249  0.374  0.499  0.624  0.749  0.875  1.000         *
'*  0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000         *
'*  0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000         *
'*  0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000         *
'*  0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000         *
'*  0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000         *
'*  0.000  0.124  0.249  0.373  0.498  0.624  0.749  0.875  1.000         *
'*  0.000  0.124  0.249  0.373  0.498  0.624  0.749  0.874  1.000         *
'*                                                                        *
'* ---------------------------------------------------------------------- *
'* REFERENCE:  "Methode de calcul numerique- 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)             *
'**************************************************************************
' Note: Here, what is new is the limit conditions involving partial derivatives dT/dy.
'       Inside domain, the classical scheme is applied to calculate the residual value at current
'       point (i,j): 4*T(i,j)-T(i-1,j)-T(i+1,j)-Ti,j-1)-T(i,j+1). Must be near zero (< ER).
'       For limit lines OA and CB, the partial derivative is approximated by:
'                       dT/dy = (T(i,j+1)-T(i,j-1))/(2*dy)
'       But the point T(i,j+1) for segment CB is out of the domain and becomes a fictitious point,
'       so we apply the following scheme:  r=4*T(i,j)-2*T(i,j-1)-T(i-1,j)-T(i+1,j).
'       For the other segment OA, it becomes: r=4*T(i,j)-2*T(i,j+1)-T(i-1,j)-T(i+1,j).
'       The exact solution for this problem is T=K*x (0<x<C), here k=0.25.
'       This allows verifying the results accuracy.
'-------------------------------------------------------------------------------------------------         
'Program laplace2

DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0

                      '           dT/dy=0
      IT = 100        '       C-------------B
      ER = .001       '       |             |
      OM = 1.3        '       |             |
      C = 4           '   T=0 |             |T=1
      NM = 8          '       |             |
                      '       |             |
                      '       O-------------A
                      '           dT/dy=0

'Labels: 10, 20

DIM T(NM + 1, NM + 1)

  'Fix limit conditions
  FOR i = 0 TO NM
    FOR j = 0 TO NM
      T(i, j) = 0#
    NEXT j
  NEXT i
  FOR j = 0 TO NM
    T(NM, j) = 1#
  NEXT j

  CLS
  n1 = NM - 1
  'main calculation loop
  FOR l = 1 TO IT
    rm = 0#
    j = 0
    FOR i = 1 TO n1
      GOSUB 1010 'call Schema1
    NEXT i
    FOR j = 1 TO n1
      FOR i = 1 TO n1
        GOSUB 1000 'call Schema
      NEXT i
    NEXT j
    j = NM
    FOR i = 1 TO n1
      GOSUB 1020 'call Schema2
    NEXT i
    IF rm < ER THEN GOTO 10
  NEXT l
  'print results
  PRINT " No convergence after "; IT; " iterations."
  PRINT " Residual rm = "; rm
  GOTO 20
10 PRINT
  PRINT " Niter  Max. Residual     W"
  PRINT USING "####"; l;
  PRINT USING "     ##.####"; rm;
  PRINT USING "       ###.##"; OM
  PRINT " Temperature"
  FOR j = NM TO 0 STEP -1
    FOR i = 0 TO NM
      PRINT USING "###.###"; T(i, j);
    NEXT i
    PRINT
  NEXT j

20 END


'calculate residual for inside square domain
1000  'Subroutine Schema
      r = T(i + 1, j) + T(i - 1, j) + T(i, j + 1) + T(i, j - 1) - 4 * T(i, j)
      T(i, j) = T(i, j) + .25 * OM * r
      IF ABS(r) > rm THEN rm = ABS(r)
RETURN

'calculate residual for lower limit
1010  'Subroutine Schema1
      r = T(i + 1, j) + T(i - 1, j) + 2 * T(i, j + 1) - 4 * T(i, j)
      T(i, j) = T(i, j) + .25 * OM * r
      IF ABS(r) > rm THEN rm = ABS(r)
RETURN

'calculate residual for upper limit
1020  'Subroutine Schema2
      r = T(i + 1, j) + T(i - 1, j) + 2 * T(i, j - 1) - 4 * T(i, j)
      T(i, j) = T(i, j) + .25 * OM * r
      IF ABS(r) > rm THEN rm = ABS(r)
RETURN

'end of file laplace2.bas