'************************************************************************ '* 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