'********************************************************************************************************* '* Solve Laplace Equation by relaxation Method: d2T/dx2 + d2T/dy2 = 0 (2) * '* ----------------------------------------------------------------------------------------------------- * '* Example #2: Determine the temperatures T(x,y) in a domain [x,y] made of two concentrical rectangles * '* (see sketch below) with following limit conditions: * '* 1) T = 0 for sides of inner rectangle. * '* 2) T = 1 for sides of outer rectangle. * '* Main parameters: * '* IM,JM: length, width of outer rectangle. * '* IB-IA,JB-JA: length, width of inner rectangle. * '* Number of subdivisions 16 in ox and 8 in Oy. * '* integration steps dx=IM/16, dy = JM/8. * '* weight coefficient w = 1.45 * '* Maximum residual error ER = 0.05 * '* Maximum number of iterations IT = 40 * '* (See file laplace.bas for example #1). * '* ----------------------------------------------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Output file laplace.lst contains: * '* * '* Niter Max. Residual W * '* 10 0.0084 1.45 * '* Temperature * '* 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 * '* * '* 1.000 0.965 0.925 0.871 0.796 0.688 0.551 0.514 0.508 0.515 0.551 0.688 0.796 0.871 0.924 0.965 1.000 * '* * '* 1.000 0.936 0.862 0.764 0.624 0.404 0.000 0.000 0.000 0.000 0.000 0.405 0.625 0.765 0.862 0.936 1.000 * '* * '* 1.000 0.917 0.822 0.700 0.533 0.304 0.000 0.000 0.305 0.535 0.700 0.822 0.917 1.000 * '* * '* 1.000 0.910 0.809 0.679 0.506 0.278 0.000 0.000 0.280 0.508 0.679 0.808 0.910 1.000 * '* * '* 1.000 0.919 0.823 0.700 0.534 0.304 0.000 0.000 0.303 0.535 0.701 0.822 0.917 1.000 * '* * '* 1.000 0.940 0.865 0.766 0.625 0.404 0.000 0.000 0.000 0.000 0.000 0.412 0.629 0.767 0.863 0.937 1.000 * '* * '* 1.000 0.965 0.928 0.872 0.797 0.688 0.551 0.514 0.507 0.514 0.550 0.688 0.803 0.874 0.926 0.966 1.000 * '* * '* 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.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 Laplace1 DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 ' Y IT = 40 ' JM --------------------------- ER = .05 ' | | OM = 1.45 ' JB|........-------- | IM = 16 ' | | T=0 | | JM = 8 ' | | | |T=1 IA = 6 ' JA|........------- | IB = 10 ' | . . | JA = 2 ' | . . | JB = 6 ' --------------------------- X ' 0 IA IB IM ' r,rm: Double DIM T(17, 9) OPEN "laplace.lst" FOR OUTPUT AS #1 'Fix limit conditions FOR i = 0 TO IM FOR j = 0 TO JM T(i, j) = 0# NEXT j NEXT i FOR i = 0 TO IM T(i, 0) = 1# T(i, JM) = 1# NEXT i FOR j = 0 TO JM T(0, j) = 1# T(IM, j) = 1# NEXT j l = 0: rm = 1# 'main calculation loop WHILE rm > ER AND l <= IT rm = 0#: l = l + 1 FOR i = 1 TO IM - 1 FOR j = 1 TO JA - 1 GOSUB 1000 'call Schema NEXT j NEXT i FOR j = JA TO JB FOR i = 1 TO IA - 1 GOSUB 1000 NEXT i FOR i = IB + 1 TO IM - 1 GOSUB 1000 NEXT i NEXT j FOR i = 1 TO IM - 1 FOR j = JB + 1 TO JM - 1 GOSUB 1000 NEXT j NEXT i WEND 'print results IF l < IT THEN PRINT #1, PRINT #1, " Niter Max. Residual W" PRINT #1, " "; l; PRINT #1, USING " ##.####"; rm; PRINT #1, USING " ####.##"; OM PRINT #1, " Temperature" FOR j = JM TO JB STEP -1 FOR i = 0 TO IM GOSUB 2000 'Print T(i,j) NEXT i PRINT #1, : PRINT #1, NEXT j FOR j = JB - 1 TO JA + 1 STEP -1 FOR i = 0 TO IA GOSUB 2000 NEXT i FOR i = IA + 1 TO IB - 1 PRINT #1, " "; NEXT i FOR i = IB TO IM GOSUB 2000 NEXT i PRINT #1, : PRINT #1, NEXT j FOR j = JA TO 0 STEP -1 FOR i = 0 TO IM GOSUB 2000 NEXT i PRINT #1, : PRINT #1, NEXT j ELSE PRINT #1, " No convergence, R max. ="; rm END IF CLOSE #1 CLS PRINT PRINT " Results in file laplace.lst..." PRINT END 'of main program 'calculate residual error at (i,j) 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 = r RETURN 'print current value of temperature array T 2000 'Subroutine Print PRINT #1, USING "##.###"; T(i, j); RETURN 'end of file laplace1.bas