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