/*********************************************************************** * 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.003131 1.450000 * * 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]. * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***********************************************************************/ #include #include // y T=x^3 #define IT 40 // ------------ #define ER 0.01 // | C | #define OM 1.45 // | | #define NM 8 // T=0| C|T=16x #define C 4 // | | // | | // ------------ x // T=0 double dx,r,rm; int i,j,l,n1; double T[NM+1][NM+1]; void Schema() { r = T[i+1][j]+T[i-1][j]+T[i][j+1]+T[i][j-1]-4.0*T[i][j]; } void main() { dx = (double) C/NM; for (i=0; i<=NM; i++) for (j=0; j<=NM; j++) T[i][j]=0.0; for (i=0; i<=NM; i++) { T[i][NM] = dx*i*dx*i*dx*i; T[NM][i] = 16.0*dx*i; } n1=NM-1; l=0; do { rm=0; l++; for (i=1; i<=n1; i++) for (j=1; j<=n1; j++) { Schema(); T[i][j] += 0.25*OM*r; if (fabs(r) > rm) rm=fabs(r); } } while (rm > ER); if (l-1; j--) { for (i=0; i<=NM; i++) printf("%8.3f",T[i][j]); printf("\n"); } } else printf("No Convergence, RMAXI=%f\n", rm); } //end of file Laplace.cpp