program Laplace c This code was written by Tom Niedzwiecki (tniedzw) to solve the c Laplace's equation on the rectangle with boundary condutions c shown below. THe Gauss-Jordan elimination code used to solve c the linear equation was taken from Numerical Recipes. c c To compile this code, name the file laplace.for, and type c c f77 -o laplace laplace.f c c The executable code will be called laplace, and its output will c appear in laplace.dat. To run the code, just type "laplace" in c the directory where the code appears. The variables used c in this code are self-explanatory, and are listed below. c c VARIABLES c c delta_x spatial incriment of grid in x-direction c delta_y spatial incriment of grid in y-direction c delta_x2 delta_x**2 c delta_y2 delta_y**2 c x_right rightmost x coordinate c x_left leftmost x coordinate c y_bottom lowermost y coordinate c y_top uppermost y coordinate c x_loc present value of x c U matrix of coefficients -- determines the spatial relat.. c V vector of potentials c Nx integer linear x-dimension c Ny integer linear y-dimension c NUM maximum number of interior unknown points c EQNUM equation number c c This code numerically solves the following two-dimensional BVP c c V2 c __________________ c | | c | | c V1 | | V3 c | | c |________________| c V4 c c c PARAMETER(NUM=400) C C PARTIAL DIFF. EQ. FOR 2D LAPLACE EQ C REAL*4 PI,delta_x,delta_x2,delta_y,delta_y2,x_right,x_left REAL*4 y_top,y_bottom REAL*4 x_loc,y_loc INTEGER ix,iy,Nx,Ny,EQNUM REAL*4 U(NUM,NUM),V(NUM) c PI=acos(-1.0) write(6,*) ' V2' write(6,*) ' __________________' write(6,*) ' | |' write(6,*) ' | |' write(6,*) ' V1 | | V3' write(6,*) ' | |' write(6,*) ' |________________|' write(6,*) ' V4' write(6,*) ' ' write(6,*) 'enter V1,V2,V3,V4' read(5,*) V1,V2,V3,V4 write(6,*) 'enter x_left, x_right, y_bottom, y_top' read(5,*) x_left, x_right, y_bottom, y_top write(6,*) 'enter Nx,Ny' read(5,*) Nx,Ny if(Nx*Ny.gt.NUM) then write(6,*) 'too fine of a grid is requested' write(6,*) 'either increase NUM and NMAX' write(6,*) 'or decrease Nx*Ny' stop end if delta_x=(x_right-x_left)/(Nx+1) delta_y=(y_top-y_bottom)/(Ny+1) c delta_x=.1 c delta_y=.1 c x_right=1. c x_left=0. c y_top=1. c y_bottom=0. C Nx=INT((x_right-x_left)/delta_x)-1 C Ny=INT((y_top-y_bottom)/delta_y)-1 delta_y2=delta_y**2 delta_x2=delta_x**2 DO I=1,NUM V(I)=0. DO J=1,NUM U(I,J)=0. ENDDO ENDDO EQNUM=0 DO iy=1,Ny DO ix=1,Nx x_loc=FLOAT(ix)*delta_x+x_left y_loc=FLOAT(iy)*delta_y+y_bottom EQNUM=EQNUM+1 c Scan through the matrix, assigning the coefficients IF(ix.EQ.1) THEN IF(iy.EQ.1) THEN c at llh corner U(EQNUM,Nx*((iy+0)-1)+(ix+1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy+1)-1)+(ix+0))= delta_x2 V(EQNUM)= -V1*delta_y2-V4*delta_x2 C ELSE IF(iy.EQ.Ny) THEN c at the ulh corner U(EQNUM,Nx*((iy+0)-1)+(ix+1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy-1)-1)+(ix+0))= delta_x2 V(EQNUM)= -V1*delta_y2-V2*delta_x2 C ELSE c at the left edge U(EQNUM,Nx*((iy+0)-1)+(ix+1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy+1)-1)+(ix+0))= delta_x2 U(EQNUM,Nx*((iy-1)-1)+(ix+0))= delta_x2 V(EQNUM)= -V1*delta_y2 C ENDIF ELSE IF (ix.EQ.Nx) THEN IF(iy.EQ.1) THEN c at the lrh corner U(EQNUM,Nx*((iy+0)-1)+(ix-1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy+1)-1)+(ix+0))= delta_x2 V(EQNUM)= -V3*(delta_x2)-V4*delta_y2 C ELSE IF(iy.EQ.Ny) THEN c at the urh corner U(EQNUM,Nx*((iy+0)-1)+(ix-1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy-1)-1)+(ix+0))= delta_x2 V(EQNUM)= -V3*(delta_x2) -V2*delta_y2 C ELSE c along the right edge U(EQNUM,Nx*((iy+0)-1)+(ix-1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy+1)-1)+(ix+0))= delta_x2 U(EQNUM,Nx*((iy-1)-1)+(ix+0))= delta_x2 V(EQNUM)= -V3*(delta_x2) C ENDIF ELSE IF(iy.EQ.1)THEN c at the lower edge U(EQNUM,Nx*((iy+0)-1)+(ix+1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix-1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy+1)-1)+(ix+0))= delta_x2 V(EQNUM)=-V4*delta_y2 C ELSE IF(iy.EQ.Ny) THEN c at the upper edge U(EQNUM,Nx*((iy+0)-1)+(ix+1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix-1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy-1)-1)+(ix+0))= delta_x2 V(EQNUM)= -V2*delta_y2 ELSE C the interior of the rectangle U(EQNUM,Nx*((iy+0)-1)+(ix+1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix-1))= delta_y2 U(EQNUM,Nx*((iy+0)-1)+(ix+0))= -2.*(delta_y2+delta_x2) U(EQNUM,Nx*((iy+1)-1)+(ix+0))= delta_x2 U(EQNUM,Nx*((iy-1)-1)+(ix+0))= delta_x2 V(EQNUM)=0.0 ENDIF ENDDO ENDDO CALL GAUSSJ(U,EQNUM,NUM,V,1,1) c Output the data WRITE(*,*) I=0 OPEN(UNIT=7,FILE='laplace.dat') DO iy=1,Ny DO ix=1,Nx I=I+1 x_loc=FLOAT(ix)*delta_x+x_left y_loc=FLOAT(iy)*delta_y+y_bottom WRITE(7,100) x_loc,y_loc,V(I) 100 FORMAT( 5X,F5.3,5X,F5.3,5X,F8.5) ENDDO ENDDO CLOSE(UNIT=7) END SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) c c Linear equation solution by Gauss-Jordan elimination. A is c an imput matrix if N by N elements, stored in an array of physical c dimensions NP by NP. B is an input matrix of N by M containing the M c right-hand side vectors, stored in an array of physical dimensions c NP by MP. On output, A is replaced by its matrix inverse, c and B is replaced by the corresponding set of solution vectors. c PARAMETER (NMAX=400) DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX) DO 11 J=1,N IPIV(J)=0 11 CONTINUE DO 22 I=1,N BIG=0. DO 13 J=1,N IF(IPIV(J).NE.1)THEN DO 12 K=1,N IF (IPIV(K).EQ.0) THEN IF (ABS(A(J,K)).GE.BIG)THEN BIG=ABS(A(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN PAUSE 'Singular matrix' ENDIF 12 CONTINUE ENDIF 13 CONTINUE IPIV(ICOL)=IPIV(ICOL)+1 IF (IROW.NE.ICOL) THEN DO 14 L=1,N DUM=A(IROW,L) A(IROW,L)=A(ICOL,L) A(ICOL,L)=DUM 14 CONTINUE DO 15 L=1,M DUM=B(IROW,L) B(IROW,L)=B(ICOL,L) B(ICOL,L)=DUM 15 CONTINUE ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.' PIVINV=1./A(ICOL,ICOL) A(ICOL,ICOL)=1. DO 16 L=1,N A(ICOL,L)=A(ICOL,L)*PIVINV 16 CONTINUE DO 17 L=1,M B(ICOL,L)=B(ICOL,L)*PIVINV 17 CONTINUE DO 21 LL=1,N IF(LL.NE.ICOL)THEN DUM=A(LL,ICOL) A(LL,ICOL)=0. DO 18 L=1,N A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 18 CONTINUE DO 19 L=1,M B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 19 CONTINUE ENDIF 21 CONTINUE 22 CONTINUE DO 24 L=N,1,-1 IF(INDXR(L).NE.INDXC(L))THEN DO 23 K=1,N DUM=A(K,INDXR(L)) A(K,INDXR(L))=A(K,INDXC(L)) A(K,INDXC(L))=DUM 23 CONTINUE ENDIF 24 CONTINUE RETURN END