!************************************************************************ !* 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]. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !************************************************************************ PROGRAM LAPLACE parameter(IT=40,ER=0.01,OM=1.45,NM=8, IC=4) ! y T=x^3 ! ------------ ! | C | ! | | ! T=0| C|T=16x ! | | ! | | ! ------------ x ! T=0 real T(0:NM,0:NM) rm=2.*ER dx = (1.0*IC)/NM T=0. do i=0, NM T(i,NM) = (dx*i)**3 T(NM,i) = 16.*dx*i end do n1=NM-1; l=0 do while (rm > ER.and.l<=IT) rm=0. l=l+1 do i=1, n1 do j=1, n1 call Schema(r,i,j,T) T(i,j) = T(i,j) + 0.25*OM*r if (abs(r) > rm) rm=abs(r) end do end do end do if (l