C ALGORITHM 621 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 4, C DEC., 1984, P. 378. C PROGRAM EXAM1(OUTPUT,TAPE6=OUTPUT) C----------------------------------------------------------------- C TEST PROGRAM 1 FOR THE ALGORITHM BDMG. C THIS TEST PROBLEM IS A TWO-DIMENSIONAL NONLINEAR PARABOLIC C EQUATION DISCUSSED IN SECTION 5 ("FIRST EXAMPLE") OF THE C COMPANION PAPER "SOFTWARE WITH LOW STORAGE REQUIREMENTS FOR C TWO-DIMENSIONAL NONLINEAR PARABOLIC DIFFERENTIAL EQUATIONS" BY C B.P. SOMMEIJER AND P.J. VAN DER HOUWEN, WHERE ALSO THE C NUMERICAL RESULTS OF THIS TEST CAN BE FOUND. C----------------------------------------------------------------- DIMENSION V(25,25),X(49),Y(49),INFO(18) DIMENSION WORK(3850),IWORK(43) EXTERNAL G,SPRAD C----------------------------------------------------------- C DEFINITION OF THE FINEST GRID AND OF THE NUMBER OF GRIDS C----------------------------------------------------------- NX = 25 NY = 25 DX = 1.0/(NX - 1.0) DO 10 I = 1,NX 10 X(I) = (I - 1)*DX DY = 1.0/(NY - 1.0) DO 20 I = 1,NY 20 Y(I) = (I - 1)*DY M = 4 C------------------------------------------------------------ C DEFINITION OF THE INTEGRATION INTERVAL AND INITIALIZATION C------------------------------------------------------------ T = 0.0 TEND = 1.0 DO 30 J = 1,NY DO 30 I = 1,NX 30 V(I,J) = SOL (T, X(I), Y(J)) C--------------------------------------- C DEFINITION OF THE PROBLEM PARAMETERS C--------------------------------------- TOL = 1.0E-4 INFO(1) = 0 INFO(2) = 2000 INFO(3) = 3 IWORK(43) = 3850 METH = 1 C-------------------------- C CALL FOR THE INTEGRATOR C-------------------------- CALL BDMG (NX, NY, M, X, Y, + T, TEND, V, G, SPRAD, TOL, METH, + WORK, IWORK, INFO, IFLAG) C---------------------------- C CHECK ERRORFLAG ON RETURN C---------------------------- IF (IFLAG .NE. 0) WRITE (6,1000) IFLAG IF (IFLAG .GE. 1 .AND. IFLAG .LE. 4) GO TO 60 C------------------------- C OUTPUT SOME STATISTICS C------------------------- WRITE (6,1001) INFO(7), INFO(8), INFO(9) WRITE (6,1002) INFO(10) C--------------------------------------------- C PRINT THE SOLUTION AND DETERMINE THE ERROR C--------------------------------------------- WRITE (6,1003) T NXM1 = NX - 1 NYM1 = NY - 1 AE = -1.0 DO 50 J = 2,NYM1 WRITE (6,1004) J DO 40 I = 2,NXM1 AE = AMAX1(AE,ABS(SOL(T, X(I), Y(J))-V(I,J))) 40 CONTINUE 50 WRITE (6,1005) (V(I,J),I = 2,NXM1) WRITE (6,1006) -ALOG10(AE) 1000 FORMAT(27H ERRORFLAG HAS BEEN SET TO:,I5) 1001 FORMAT(27H NUMBER OF STEPS PERFORMED:,I6/ + 26H NUMBER OF REJECTED STEPS:,I6/ + 31H NUMBER OF STEPS WITH OVERFLOW:,I6) 1002 FORMAT(44H TOTAL NUMBER OF (FINEST GRID)G-EVALUATIONS:,I6) 1003 FORMAT(//15H SOLUTION AT T=,F10.4/) 1004 FORMAT(I4) 1005 FORMAT(5(E20.10)) 1006 FORMAT(34H MINIMAL NUMBER OF CORRECT DIGITS:,F7.2) C------------------------------------ C END OF MAIN PROGRAM FOR EXAMPLE 1 C------------------------------------ 60 STOP END FUNCTION SOL (T, X, Y) C----------------------------------------------------- C ANALYTIC FUNCTION USED FOR INITIALIZATION,BOUNDARY C CONDITIONS AND CALCULATION OF THE ERROR C----------------------------------------------------- SOL = (0.8*(2.0*T + X+ Y ))**0.25 RETURN END FUNCTION SPRAD (K, NXK, NYK, XK, YK, T, VK, TAU) DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK) C--------------------------------------------------------- C DEFINE AN UPPER ESTIMATE OF THE SPECTRAL RADIUS ON THE C INTERVAL [T,T+TAU] C--------------------------------------------------------- SPRAD = 32.0*(1.0 + T + TAU)*((NXK - 1)**2 + (NYK - 1)**2) RETURN END SUBROUTINE G (K, NXK, NYK, XK, YK, T, VK, DVK) DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK),DVK(NXK,NYK) C---------------------------------------------- C DEFINE THE DERIVATIVE FUNCTION G AT LEVEL K C---------------------------------------------- NXKM1 = NXK - 1 NYKM1 = NYK - 1 C----------------------------------- C TREATMENT OF THE BOUNDARY POINTS C----------------------------------- DO 10 J = 1,NYK DO 10 I = 1,NXK,NXKM1 DVK(I,J) = 0.0 10 VK(I,J) = SOL (T, XK(I), YK(J)) DO 20 J = 1,NYK,NYKM1 DO 20 I = 2,NXKM1 DVK(I,J) = 0.0 20 VK(I,J) = SOL (T, XK(I), YK(J)) C--------------------------------- C DERIVATIVES AT INTERNAL POINTS C--------------------------------- DO 30 J = 2,NYKM1 DO 30 I = 2,NXKM1 30 DVK(I,J)=(VK(I+1,J)**5-2.0*VK(I,J)**5+VK(I-1,J)**5)*NXKM1**2 + + (VK(I,J+1)**5-2.0*VK(I,J)**5+VK(I,J-1)**5)*NYKM1**2 RETURN END C PROGRAM EXAM2 (OUTPUT,TAPE6=OUTPUT) C----------------------------------------------------------------- C TEST PROGRAM 2 FOR THE ALGORITHM BDMG. C THIS TEST PROBLEM IS A TWO-DIMENSIONAL NONLINEAR PARABOLIC C EQUATION INCLUDING A MIXED-DERIVATIVE TERM, DISCUSSED IN C IN SECTION 5 ("SECOND EXAMPLE") OF THE COMPANION PAPER C "SOFTWARE WITH LOW STORAGE REQUIREMENTS FOR TWO-DIMENSIONAL C NONLINEAR PARABOLIC DIFFERENTIAL EQUATIONS" BY B.P. SOMMEIJER C AND P.J. VAN DER HOUWEN, WHERE ALSO THE NUMERICAL RESULTS OF C THIS TEST CAN BE FOUND. C----------------------------------------------------------------- DIMENSION V(17,17),X(31),Y(31),INFO(16) DIMENSION WORK(1776),IWORK(33) EXTERNAL G,SPRAD C----------------------------------------------------------- C DEFINITION OF THE FINEST GRID AND OF THE NUMBER OF GRIDS C----------------------------------------------------------- NX = 17 NY = 17 DX = 1.0/(NX - 1.0) DO 10 I = 1,NX 10 X(I) = (I - 1)*DX DY = 1.0/(NY - 1.0) DO 20 I = 1,NY 20 Y(I) = (I - 1)*DY M = 3 C--------------------------------------------------------------- C SUCCESSIVE INTEGRATION PROCESSES FOR DIFFERENT VALUES OF TOL C--------------------------------------------------------------- DO 60 ITOL = 1,6 C------------------------------------------------------------ C DEFINITION OF THE INTEGRATION INTERVAL AND INITIALIZATION C------------------------------------------------------------ T = 0.0 TEND = 1.0 DO 30 J = 1,NY DO 30 I = 1,NX 30 V(I,J) = SOL (T, X(I), Y(J)) C--------------------------------------- C DEFINITION OF THE PROBLEM PARAMETERS C--------------------------------------- TOL = 10.0**(-ITOL) INFO(1) = 0 INFO(2) = 1000 INFO(3) = 3 IWORK(33) = 1776 METH = 1 C-------------------------- C CALL FOR THE INTEGRATOR C-------------------------- CALL BDMG (NX, NY, M, X, Y, + T, TEND, V, G, SPRAD, TOL, METH, + WORK, IWORK, INFO, IFLAG) C---------------------------- C CHECK ERRORFLAG ON RETURN C---------------------------- IF (IFLAG .NE. 0) WRITE (6,1000) IFLAG IF (IFLAG .GE. 1 .AND. IFLAG .LE. 4) GO TO 60 C------------------------- C OUTPUT SOME STATISTICS C------------------------- WRITE (6,1007) TOL WRITE (6,1001) INFO(7), INFO(8), INFO(9) WRITE (6,1002) INFO(10) C--------------------------------------------- C PRINT THE SOLUTION AND DETERMINE THE ERROR C--------------------------------------------- WRITE (6,1003) T NXM1 = NX - 1 NYM1 = NY - 1 AE = -1.0 DO 50 J = 2,NYM1 WRITE (6,1004) J DO 40 I = 2,NXM1 AE = AMAX1(AE,ABS(SOL (T, X(I), Y(J))-V(I,J))) 40 CONTINUE 50 WRITE (6,1005) (V(I,J),I = 2,NXM1) WRITE (6,1006) -ALOG10(AE) 60 CONTINUE 1000 FORMAT(27H ERRORFLAG HAS BEEN SET TO:,I5) 1001 FORMAT(27H NUMBER OF STEPS PERFORMED:,I6/ + 26H NUMBER OF REJECTED STEPS:,I6/ + 31H NUMBER OF STEPS WITH OVERFLOW:,I6) 1002 FORMAT(44H TOTAL NUMBER OF (FINEST GRID)G-EVALUATIONS:,I6) 1003 FORMAT(//15H SOLUTION AT T=,F10.4/) 1004 FORMAT(I4) 1005 FORMAT(5(E20.10)) 1006 FORMAT(34H MINIMAL NUMBER OF CORRECT DIGITS:,F7.2) 1007 FORMAT(1H1,7H TOL = ,E10.2) C------------------------------------ C END OF MAIN PROGRAM FOR EXAMPLE 2 C------------------------------------ STOP END FUNCTION SOL (T, X, Y) C----------------------------------------------------- C ANALYTIC FUNCTION USED FOR INITIALIZATION,BOUNDARY C CONDITIONS AND CALCULATION OF THE ERROR C----------------------------------------------------- SOL = EXP(-T)*X*Y*(X + Y) RETURN END FUNCTION SPRAD (K, NXK, NYK, XK, YK, T, VK, TAU) DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK) C--------------------------------------------------------- C DEFINE AN UPPER ESTIMATE OF THE SPECTRAL RADIUS ON THE C INTERVAL [T,T+TAU] C--------------------------------------------------------- SPRAD = 6.0*((NXK-1)**2 + (NYK-1)**2) + 2.0*(NXK-1)*(NYK-1) RETURN END SUBROUTINE G (K, NXK, NYK, XK, YK, T, VK, DVK) DIMENSION XK(NXK),YK(NYK),VK(NXK,NYK),DVK(NXK,NYK) C---------------------------------------------- C DEFINE THE DERIVATIVE FUNCTION G AT LEVEL K C---------------------------------------------- NXKM1 = NXK - 1 NYKM1 = NYK - 1 C----------------------------------- C TREATMENT OF THE BOUNDARY POINTS C----------------------------------- DO 10 J = 1,NYK DO 10 I = 1,NXK,NXKM1 DVK(I,J) = 0.0 10 VK(I,J) = SOL (T, XK(I), YK(J)) DO 20 J = 1,NYK,NYKM1 DO 20 I = 2,NXKM1 DVK(I,J) = 0.0 20 VK(I,J) = SOL (T, XK(I), YK(J)) C--------------------------------- C DERIVATIVES AT INTERNAL POINTS C--------------------------------- ET = EXP(-T) DO 30 J = 2,NYKM1 Y2 = YK(J)**2 DO 30 I = 2,NXKM1 X2 = XK(I)**2 A = 0.5*X2 + Y2 B = -0.5*(X2 + Y2) C = X2 + 0.5*Y2 D = 1.0/(1.0 + (X2*YK(J) + XK(I)*Y2)*ET) DVK(I,J) = A*(VK(I-1,J)-2.0*VK(I,J)+VK(I+1,J))*NXKM1*NXKM1 + + C*(VK(I,J-1)-2.0*VK(I,J)+VK(I,J+1))*NYKM1*NYKM1 + + (VK(I+1,J+1)-VK(I-1,J+1)+VK(I-1,J-1)-VK(I+1,J-1))* + 2.0*B*NXKM1*NYKM1/4.0 30 DVK(I,J) = DVK(I,J)*(D*(1.0 + VK(I,J)))**10 RETURN END SUBROUTINE BDMG (NX, NY, M, X, Y, * T, TEND, V, G, SPRAD, TOL, METH, * WORK, IWORK, INFO, IFLAG) C C----------------------------------------------------------------------- C C PURPOSE AND METHOD C C----------------------------------------------------------------------- C C BDMG IS A TIME-INTEGRATOR DESIGNED TO SOLVE A SYSTEM OF ORDINARY C DIFFERENTIAL EQUATIONS(ODE'S) ORIGINATING FROM SEMI-DISCRETIZATION C OF A (SCALAR) GENERAL PARABOLIC PARTIAL DIFFERENTIAL EQUATION(PDE) C IN TWO SPACE DIMENSIONS. THIS SEMI-DISCRETIZATION CAN BE PERFORMED C EITHER BY HAND OR BY PDETWO[3], AN INTERFACE TO FORM AND EVALUATE C A SEMI-DISCRETE APPROXIMATION TO THE ORIGINAL PDE. C C THE METHOD ON WHICH BDMG IS BASED [2], CONSISTS OF THE APPLICATION C OF THE SECOND ORDER BACKWARD DIFFERENTIATION FORMULA TO THIS SEMI- C DISCRETE SYSTEM. THE RESULTING NON-LINEAR PROBLEM IS SOLVED USING A C MULTIGRID TECHNIQUE. C C THE MAIN CHARACTERISTIC OF THE CODE BDMG IS ITS MINIMAL STORAGE C REQUIREMENTS. C C THE APPLICABILITY OF BDMG IS RESTRICTED TO C (1) TWO SPACE DIMENSIONS C (2) RECTANGULAR DOMAINS C (3) SCALAR PARTIAL DIFFERENTIAL EQUATIONS C (4) PDE'S IN WHICH DIFFUSION GREATLY DOMINATES CONVECTION C (5) PDE'S WITHOUT MIXED DERIVATIVES(THIS RESTRICTION ONLY APPLIES C IF PDETWO IS USED) C C----------------------------------------------------------------------- C C PROBLEM DEFINITION C C----------------------------------------------------------------------- C C BDMG APPLIES TO SEMI-DISCRETE SYSTEMS OF THE FORM C C C DV C (1) [ -- = G (T, V) ] ,I=1,...,NX C DT I,J J=1,...,NY , C C C WHERE THE (SCALAR) MATRIX ELEMENTS V(I,J) ARE ASSOCIATED TO THE C GRID POINTS (X(I),Y(J)) WHICH FORM A GRID WITH RECTANGULAR MESHES C IN THE (X,Y)-PLANE, INCLUDING THE BOUNDARY POINTS. C THIS SYSTEM IS SUPPOSED TO ORIGINATE FROM PARTIAL DIFFERENTIAL C EQUATIONS OF THE FORM C C C (2) UT = F (T, X, Y, U, UX, UY, UXX, UYY, UXY) , C C WHERE UT IS THE FIRST PARTIAL DERIVATIVE OF U WITH RESPECT TO T C AND UX,UXX ARE THE FIRST AND SECOND PARTIAL DERIVATIVE OF U WITH C RESPECT TO X C ETC. C THUS, ANY PDE OF THE FORM (2) DEFINED ON A RECTANGLE WITH INITIAL C CONDITION DEFINED ON THIS RECTANGLE AND BOUNDARY CONDITIONS OF C THE FORM C C (3) A(T,X,Y)*U + B(T,X,Y)*UN = C(T,X,Y), UN IS DERIVATIVE OF U C NORMAL TO THE BOUNDARY C C CAN BE DEALT WITH BY BDMG AS SOON AS (2) AND (1) ARE SEMI- C DISCRETIZED ON THE GRID (X(I),Y(J)). C C THE INITIAL CONDITION BEING DEFINED FOR EACH V(I,J) NEEDS NOT BE C CONSISTENT WITH THE BOUNDARY CONDITIONS. C C----------------------------------------------------------------------- C C USER-REQUIRED INFORMATION C C----------------------------------------------------------------------- C C IN DEFINING HIS PROBLEM, THE USER IS ASKED TO PROVIDE A MAIN C PROGRAM IN WHICH THE GRID IS DEFINED, THE INTEGRATION INTERVAL IS C SPECIFIED AND THE INITIAL CONDITION IS SET. BECAUSE BDMG IS A C TIME-INTEGRATOR FOR A SYSTEM OF ODE'S IT REQUIRES IN ADDITION A C SEMI-DISCRETE APPROXIMATION TO THE ORIGINAL PDE. THIS SEMI- C DISCRETIZATION CAN BE PERFORMED EITHER BY HAND OR BY PDETWO; C BY MEANS OF THE PARAMETER METH THE USER CAN INDICATE HIS CHOICE. C SELECTING THE FIRST POSSIBILITY, A SUBROUTINE IS REQUIRED C DELIVERING THE DERIVATIVE FUNCTION G IN (1). IN CASE OF THE USAGE C OF PDETWO IT SUFFICES TO SUPPLY ROUTINES DEFINING THE PARTIAL C DIFFERENTIAL EQUATION AND THE BOUNDARY CONDITIONS (SEE [3] FOR A C DETAILED DESCRIPTION OF THE USER-REQUIRED INFORMATION). C C----------------------------------------------------------------------- C C THE PARAMETERS C C----------------------------------------------------------------------- C C WE DISTINGUISH BETWEEN THREE TYPES OF PARAMETERS: C C GRID-DEFINING PARAMETERS C ------------------------ C C NX,NY ARE THE NUMBER OF MESH LINES IN X- AND Y-DIRECTION, C RESPECTIVELY C M IS THE NUMBER OF SUCCESSIVE GRIDS, USED IN THE C MULTIGRID METHOD(FOR A DISCUSSION OF M SEE BELOW). C X(1),...,X(NX) CONTAIN THE X-COORDINATES OF THE MESH POINTS C ALONG EACH GRID LINE IN THE X-DIRECTION. C (X(1) .LT. X(2) .LT. ... .LT. X(NX)). C Y(1),...,Y(NY) CONTAIN THE ANALOGOUS INFORMATION IN THE C Y-DIRECTION. C (Y(1) .LT. Y(2) .LT. ... .LT. Y(NY)). C C NX,NY,M,X,Y ARE INPUT PARAMETERS AND HAVE TO BE INITIALIZED BY C THE USER. C NOTE THAT X AND Y ARE NOT NECESSARILY EQUIDISTANT. C C IT IS WORTH GOING INTO MORE DETAIL ABOUT THESE PARAMETERS. C BECAUSE THE METHOD IS BASED ON A MULTIGRID METHOD, A SEQUENCE OF M C SUCCESSIVE GRIDS IS NECESSARY; THESE GRIDS ARE GENERATED BY THE C CODE. THE LEVEL INDEX K RUNS FROM 1 (WHICH IS ASSOCIATED WITH THE C COARSEST GRID) UNTIL M (THE FINEST GRID, WHICH IS IDENTICAL TO THE C GRID CHOOSEN BY THE USER AND DEFINED BY MEANS OF THE PARAMETERS C NX,NY,X,Y). INTERNAL VALUES NXK, NYK ARE CALCULATED AS WELL AS C ARRAYS XK(I),I=1,...,NXK AND YK(J),J=1,...,NYK, DEFINING THE GRID C AT LEVEL K. THESE VALUES ARE USED, AMONG OTHER THINGS, TO PASS TO C THE SUBROUTINE G (SEE BELOW) IN WHICH THE USER HAS TO DEFINE THE C DERIVATIVE AT THAT PARTICULAR GRID (XK(I),YK(J)) (ONLY IF METH=1, C SEE BELOW). C C A GRID AT LEVEL K IS OBTAINED FROM A GRID AT LEVEL K+1 (DENOTED C BY KP1) BY ALTERNATINGLY DELETING GRID LINES;HENCE, NXK=(NXKP1+1)/2 C XK(1)=XKP1(1), XK(2)=XKP1(3), ETC. THIS MEANS THAT THE GRID LINES C ON LEVEL K COINCIDE WITH THOSE ON LEVEL K+1. BECAUSE THE VALUES OF C THE ARRAYS XK (K