/******************************************************************************************************** * Solve Laplace Equation by relaxation Method: d2T/dx2 + d2T/dy2 = 0 (2) * * ----------------------------------------------------------------------------------------------------- * * Example #2: Determine the temperatures T(x,y) in a domain [x,y] made of two concentrical rectangles * * (see sketch below) with following limit conditions: * * 1) T = 0 for sides of inner rectangle. * * 2) T = 1 for sides of outer rectangle. * * Main parameters: * * IM,JM: length, width of outer rectangle. * * IB-IA,JB-JA: length, width of inner rectangle. * * Number of subdivisions 16 in ox and 8 in Oy. * * integration steps dx=IM/16, dy = JM/8. * * weight coefficient w = 1.45 * * Maximum residual error ER = 0.05 * * Maximum number of iterations IT = 40 * * (See file laplace.pas for example #1). * * ----------------------------------------------------------------------------------------------------- * * SAMPLE RUN: * * * * Output file laplace.lst contains: * * * * Niter Max. Residual W * * 10 0.008388 1.45 * * Temperature * * 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 * * * * 1.000 0.965 0.925 0.871 0.796 0.688 0.551 0.514 0.508 0.515 0.551 0.688 0.796 0.871 0.924 0.965 1.000 * * * * 1.000 0.936 0.862 0.764 0.624 0.404 0.000 0.000 0.000 0.000 0.000 0.405 0.625 0.765 0.862 0.936 1.000 * * * * 1.000 0.917 0.822 0.700 0.533 0.304 0.000 0.000 0.305 0.535 0.700 0.822 0.917 1.000 * * * * 1.000 0.910 0.809 0.679 0.506 0.278 0.000 0.000 0.280 0.508 0.679 0.808 0.910 1.000 * * * * 1.000 0.919 0.823 0.700 0.534 0.304 0.000 0.000 0.303 0.535 0.701 0.822 0.917 1.000 * * * * 1.000 0.940 0.865 0.766 0.625 0.404 0.000 0.000 0.000 0.000 0.000 0.412 0.629 0.767 0.863 0.937 1.000 * * * * 1.000 0.965 0.928 0.872 0.797 0.688 0.551 0.514 0.507 0.514 0.550 0.688 0.803 0.874 0.926 0.966 1.000 * * * * 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 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) * ********************************************************************************************************/ #include #include const // Y IT=40, // JM --------------------------- IM=16, // | | JM=8, // JB|........-------- | IA=6, // | | T=0 | | IB=10, // | | | |T=1 JA=2, // JA|........------- | JB=6; // | . . | // | . . | // --------------------------- X // 0 IA IB IM double r,rm; int i,j,k,l; double T[17][9]; double ER = 0.05; double OM = 1.45; FILE *fp; //calculate residual error at (i,j) 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=r; } //print current value of temperature array T void Print() { fprintf(fp,"%6.3f", T[i][j]); } void main() { fp = fopen("laplace.lst","w"); //Fix limit conditions for (i=0; i<=IM; i++) for (j=0; j<=JM; j++) T[i][j]=0.0; for (i=0; i<=IM; i++) { T[i][0]=1.0; T[i][JM]=1.0; } for (j=0; j<=JM; j++) { T[0][j]=1.0; T[IM][j]=1.0; } l=0; //main calculation loop do { rm=0.0; l++; for (i=1; i=ER && l <= IT); //print results if (l < IT) { fprintf(fp,"\n Niter Max. Residual W\n"); fprintf(fp,"%4d %f %6.2f\n", l, rm, OM); fprintf(fp," Temperature\n"); for (j=JM; j>=JB; j--) { for (i=0; i<=IM; i++) Print(); fprintf(fp,"\n"); fprintf(fp,"\n"); } for (j=JB-1; j>=JA+1; j--) { for (i=0; i<=IA; i++) Print(); for (i=IA+1; i<=IB-1; i++) fprintf(fp," "); for (i=IB; i<=IM; i++) Print(); fprintf(fp,"\n"); fprintf(fp,"\n"); } for (j=JA; j>=0; j--) { for (i=0; i<=IM; i++) Print(); fprintf(fp,"\n"); fprintf(fp,"\n"); } } else fprintf(fp," No convergence, R max. = %f\n", rm); fclose(fp); printf("\n See results in file laplace.lst...\n\n"); } // end of file laplace1.cpp