/************************************************************************* * 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.000860 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: "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) * ************************************************************************** { 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 #include const // dT/dy=0 IT = 100, // C-------------B C = 4, // | | NM = 8; // | | // T=0 | |T=1 // | | // | | // O-------------A // dT/dy=0 //Label: e10 double T[NM+1][NM+1]; int i,j,l,n1; double r, rm=1.0; double ER=0.001, OM=1.3; //calculate residual for inside square domain void 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]+0.25*OM*r; if (fabs(r) > rm) rm=fabs(r); } //calculate residual for lower limit void 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]+0.25*OM*r; if (fabs(r) > rm) rm=fabs(r); } //calculate residual for upper limit void 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]+0.25*OM*r; if (fabs(r) > rm) rm=fabs(r); } void main() { //Fix limit conditions for (i=0; i<=NM; i++) for (j=0; j<=NM; j++) T[i][j]=0.0; for (j=0; j<=NM; j++) T[NM][j]=1.0; n1=NM-1; //main calculation loop for (l=1; l<=IT; l++) { rm=0.0; j=0; for (i=1; i<=n1; i++) Schema1(); for (j=1; j<=n1; j++) for (i=1; i<=n1; i++) Schema(); j=NM; for (i=1; i<=n1; i++) Schema2(); if (rm < ER) goto e10; } //print results printf("\n No convergence after %d iterations.\n", IT); printf(" Residual rm = %f\n", rm); return; e10:; printf("\n Niter Max. Residual W\n"); printf("%4d %f %6.2f\n", l, rm, OM); printf(" Temperature\n"); for (j=NM; j>=0; j--) { for (i=0; i<=NM; i++) printf("%7.3f", T[i][j]); printf("\n"); } printf("\n"); } // end of file laplace2.cpp