C ALGORITHM 685, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 16, NO. 4, PP. 325-351. C-----TOMS SERRG2 (MAIN PORTION OF ALGORITHM) SUBROUTINE SERRG2(K, M, X, N, Y, P1, R, Q1, P2, Q2, F, V, 1 IA, BETAA, GAMA, IB, BETAB, GAMB, IC, BETAC, GAMC, ID, BETAD, 1 GAMD) C C SERRG2 - SEPARABLE ELLIPTIC RAYLEIGH-RITZ-GALERKIN TWO C DIMENSIONAL PARTIAL DIFFERENTIAL EQUATION SOLVER. C C THIS SUBROUTINE IS DESIGNED TO SOLVE THE SEPARABLE ELLIPTICAL C P.D.E.PROBLEM C (L(X) + M(Y))U(X,Y)=F(X,Y) (1) C WHERE C L(X)= -D (P1(X)D ) + R(X)D + Q1(X) C X X X C C AND M(Y) = -D (P2(Y)D )+Q2(Y) C Y Y C C ON THE RECTANGLE A<=X<=B, C<=Y<=D. C C THE SUBROUTINE RETURNS THE COEFFICIENTS V TO A RAYLEIGH-RITZ- C GALERKIN APPROXIMATION TO U OF THE FORM C C M N C U(X,Y)=SUM ( SUM V(I,J) B (X) B(Y)) (2) C I=1 J=1 I J C C WHERE THE B(X)'S ARE B-SPLINES OF ORDER K ASSOCIATED WITH THE MESH C A=X =.=X < X <=.....<=X 2 C 120 CONTINUE VALKNT = 0.5D0 CALL BKMSH(C, D, MESH, K, VALKNT, Y, N) CALL SERRG2(K, M, X, N, Y, ONE, ZERO, ZERO, P212, ZERO, 1 F12, V, 1, ZR,GAMA12,1,ZR,GAMB12,1,ZR,GAMC12,1,ZR,GAMD12) CALL EREST3(K, M, X,0,N, Y,0,V, UTRU12, ERROR) GOTO 1000 C C ALL DIRICHLET - MULTIPLE KNOTS IN Y FOR K > 2 C 130 CONTINUE VALKNT = 0.5D0 CALL BKMSH(C, D, MESH, K, VALKNT, X, M) CALL SERRG2(K, M, X, N, Y, P212, ZERO, ZERO, ONE, ZERO, 1 F13, V, 1, ZR,GAMC12,1,ZR,GAMD12,1,ZR,GAMA12,1,ZR,GAMB12) CALL EREST3(K, M, X,0,N, Y,0,V, UTRU13, ERROR) GOTO 1000 C C ALL DIRICHLET - MULTIPLE KNOTS IN X AND Y FOR K > 2 C 140 IA = 1 IB = 1 IC = 1 ID = 1 VALKNT = 0.5D0 CALL BKMSH(C, D, MESH, K, VALKNT, Y, N) CALL BKMSH(C, D, MESH, K, VALKNT, X, M) C CALL SERRG2(K, M, X, N, Y, P212, ZERO, ZERO, P212, ZERO, 1 F14, V, IA, ZR,GAMA12,IB,ZR,GAMB14,IC,ZR,GAMA12,ID,ZR,GAMB14) CALL EREST3(K, M, X,0,N, Y,0,V, UTRU14, ERROR) GOTO 1000 C 1000 WRITE(OUTUNT,51)IPROB,K,MESH,ERROR(1),ERROR(2),ERROR(3) 51 FORMAT(1H ,3I9,3E12.3) GO TO 10 C 2000 STOP END SUBROUTINE UMESH(A, B, NX, K, T, N) C THIS SUBROUTINE COMPUTES A UNIFORM MESH WITH NX-1 INTERVALS C WITH THE ENDPOINTS AT A AND B HAVING MULTIPLICITY K C C INPUT PARAMETERS C C A THE VALUE OF THE FIRST POINT IN THE MESH C B THE VALUE OF THE LAST POINT OF THE MESH C NX THE NUMBER OF INTERVALS -1 IN THE UNIFORM MESH C K NUMBER OF MULTIPLE KNOTS AT THE END POINTS C C OUTPUT PARAMETERS C C T THE KNOT SEQUENCE C N K+NX-2, THE NUMBER OF B-SPLINES INTEGER NX, K, N DOUBLE PRECISION A, B, T(1), DX C C LEFT HAND ENDPOINTS DO 10 I=1,K 10 T(I) = A DX = (B - A)/DBLE(FLOAT(NX-1)) C COMPUTE THE INTERIOR KNOTS DO 20 I=1,NX IK1 = I + K - 1 20 T(IK1) = DBLE(FLOAT(I-1))*DX + A N = K + NX - 2 C RIGHT HAND ENDPOINTS DO 30 I=1,K NPI = N+I 30 T(NPI) = B RETURN END SUBROUTINE BKMSH(A, B, NX, K, BRK, T, N) C THIS SUBROUTINE CREATES A MESH WHICH HAS ENDPOINTS C A AND B AND A KNOT OF MULTIPLICITY K-1 AT POINT BRK C AND HAS UNIFORM INTERVALS BEFORE AND AFTER BRK WITH C NX KNOTS SPREAD PROPROTIONATELY BEFORE AND AFTER BRK C C INPUT PARAMETERS C C A THE VALUE OF THE FIRST POINT IN THE MESH C B THE VALUE OF THE LAST POINT OF THE MESH C NX THE NUMBER OF POINTS IN THE UNIFORM MESH C K NUMBER OF MULTIPLE KNOTS AT THE END POINTS C BRK THE POINT OF THE MULTIPLE KNOT C OUTPUT PARAMETERS C T THE KNOT SEQUENCE C N T HAS N+K ENTRIES INTEGER NX, K, N DOUBLE PRECISION A, B, T(1), DX DOUBLE PRECISION BRK, DX2 C C LEFT HAND ENDPOINTS DO 10 I=1,K 10 T(I) = A C NX1=(BRK-A)/(B-A)*DBLE(FLOAT(NX))+1.D0 NX2 = NX-NX1+1 IF (NX1.LT.1) GO TO 30 C C DO INTERIOR POINTS IN THE FIRST INTERVAL C DX=(BRK-A)/DBLE(FLOAT(NX1-1)) DO 20 I=1,NX1 IK=I+K-1 T(IK)=DBLE(FLOAT(I-1))*DX+A 20 CONTINUE C C DO THE BREAKPOINT C 30 N=NX1+K-2 KM1= K-1 DO 40 I=1,KM1 NPI=N+I 40 T(NPI)=BRK C C DO THE POINTS IN THE SECOND INTERVAL C IF (NX2.LT.1) GO TO 60 DX2= (B-BRK)/DBLE(FLOAT(NX2-1)) NMK2=N+K-2 DO 50 I=1,NX2 IK2 = NMK2+I T(IK2)=BRK+DBLE(FLOAT (I-1))*DX2 50 CONTINUE C C DO THE RIGHT HAND ENDPOINT 60 N=NMK2+NX2-1 DO 70 I=1,K NPI=N+I T(NPI)=B 70 CONTINUE RETURN END DOUBLE PRECISION FUNCTION FOUR(X) C THIS SUBROUTINE RETURNS THE VALUE 4 FOR EVERY X DOUBLE PRECISION X FOUR=4.0D0 RETURN END DOUBLE PRECISION FUNCTION ONE(X) C THIS SUBROUTINE RETURNS THE VALUE 1 FOR EVERY X DOUBLE PRECISION X C ONE = 1.0D0 RETURN END DOUBLE PRECISION FUNCTION ZERO(X) C THIS SUBROUTINE RETURNS THE VALUE 0 FOR EVERY X DOUBLE PRECISION X C ZERO = 0.0D0 RETURN END SUBROUTINE F5(X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR PROBLEM 1,2,3, AND 4 C DOUBLE PRECISION X, Y(LCK), FXY(LCK), T1, T2, T3 DOUBLE PRECISION SAVE(1000) INTEGER INIT COMMON /INSAV/ SAVE, INIT IF (INIT.GT.0)GO TO 3 DO 2 I=1,LCK 2 SAVE(I)=DEXP(Y(I))*Y(I) INIT=1 3 CONTINUE C T1 = -3.0D0 * DEXP(X)*X T2 = X + 3.0D0 T3 = X - 1.0D0 DO 1 I=1,LCK 1 FXY(I) = T1 * SAVE(I) * ((Y(I) - 1.0D0) * T2 + 1 T3 * (Y(I) + 3.0D0)) RETURN END C SUBROUTINE UTRUE5(X, Y,UTRU) C C IN UTRU(1) THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR PROBLEMS 1,2,3, AND 4 C IN UTRUE(2) ITS DERIVATIVE WITH RESPECT TO X IS C RETURNED AND IN UTRUE(3) ITS DERIVATIVE WITH RESPECT C TO Y IS RETURNED C DOUBLE PRECISION X, Y,UTRU(3) C UTRU(1)= 3.0D0 * DEXP(X + Y) * X * (X - 1.0D0) * Y * (Y - 1.0D0) C UTRU(2)=3.D0*DEXP(X+Y)*Y*(Y-1.D0)*(2.D0*X-1.0)+UTRU(1) UTRU(3)=3.D0*DEXP(X+Y)*X*(X-1.0D0)*(2.D0*Y-1.D0)+UTRU(1) RETURN END SUBROUTINE F9(X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR PROBLEMS 5,6, AND 8 C DOUBLE PRECISION X, Y(LCK), FXY(LCK) C DOUBLE PRECISION TWOPI, PI2, X2, C2PY(1000),PI2P1,PIO4 INTEGER INIT COMMON /INSAV/ C2PY, INIT DATA PIO4 /0.0D0/ IF (PIO4.EQ.0.0D0) PIO4=DATAN(1.0D0) PI2 = 16.D0*PIO4**2 PI2P1=PI2+1 IF (INIT.GT.0) GO TO 3 TWOPI=8.D0*PIO4 DO 2 I=1,LCK 2 C2PY(I)=DCOS(TWOPI*Y(I)) INIT=1 3 CONTINUE C X2 = 2.D0-4.D0*X*X*PI2P1 DO 1 I=1,LCK 1 FXY(I) = -C2PY(I)*X2 C RETURN END C SUBROUTINE UTRUE9(X, Y,UTRU) C C IN UTRU(1) THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR PROBLEMS 5,6, and 8 C IN UTRUE(2) ITS DERIVATIVE WITH RESPECT TO X IS C RETURNED AND IN UTRUE(3) ITS DERIVATIVE WITH RESPECT C TO Y IS RETURNED C DOUBLE PRECISION X, Y,UTRU(3) C DOUBLE PRECISION TWOPI DATA TWOPI/0.0D0/ IF (TWOPI .EQ. 0.0D0) TWOPI = 8.D0*DATAN(1.0D0) UTRU(1) = X*X*(DCOS(TWOPI*Y)) UTRU(2)=2.D0*X*DCOS(TWOPI*Y) UTRU(3)=-X*X*DSIN(TWOPI*Y)*TWOPI C RETURN END SUBROUTINE F11(X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR PROBLEMS 7 AND 9 C DOUBLE PRECISION X, Y(LCK), FXY(LCK) C DOUBLE PRECISION TWOPI, PI2, Y2, C2PX, PIO4, PI2P1 C DATA PIO4 /0.0D0/ IF (PIO4.EQ.0.0D0) PIO4=DATAN(1.0D0) PI2 = 16.D0*PIO4**2 PI2P1=PI2+1 TWOPI=8.D0*PIO4 C2PX=DCOS(TWOPI*X) DO 1 I=1,LCK Y2=Y(I)*Y(I) 1 FXY(I) = -C2PX*(2.d0-4.d0*Y2*(PI2P1)) C RETURN END C SUBROUTINE UTRU11(Y, X,UTRU) C C IN UTRU(1) THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR PROBLEMS 7 AND 9 C IN UTRUE(2) ITS DERIVATIVE WITH RESPECT TO X IS C RETURNED AND IN UTRUE(3) ITS DERIVATIVE WITH RESPECT C TO Y IS RETURNED C DOUBLE PRECISION X, Y,UTRU(3) C DOUBLE PRECISION TWOPI DATA TWOPI/0.0D0/ IF (TWOPI .EQ. 0.D0) TWOPI = 8.D0*DATAN(1.0D0) UTRU(1) = X*X*(DCOS(TWOPI*Y)) UTRU(3) = 2.D0*X*DCOS(TWOPI*Y) UTRU(2) = -X*X*TWOPI*DSIN(TWOPI*Y) C RETURN END SUBROUTINE F12 (X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR PROBLEM 10 C DOUBLE PRECISION X, Y(LCK), FXY(LCK), C2PY(1000), CONST INTEGER INIT COMMON /INSAV/ C2PY, INIT C IF (INIT .GT. 0) GO TO 3 CONST = 1.0D0 / DEXP(0.5D0) DO 2 I = 1, LCK IF (Y(I) .LE. 0.5D0) GO TO 1 C2PY(I) = - CONST * DEXP (Y(I)) GO TO 2 1 C2PY(I) = 0.0D0 2 CONTINUE INIT = 1 3 CONTINUE C DO 6 I = 1, LCK 6 FXY(I) = C2PY(I) C RETURN END SUBROUTINE F13 (X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR PROBLEM 11 C DOUBLE PRECISION X, Y(LCK), FXY(LCK), CONST C IF (X .LE. 0.5D0) GO TO 1 CONST = - 1.0D0 / DEXP(0.5D0) * DEXP(X) GO TO 2 1 CONST = 0.0D0 2 CONTINUE C DO 3 I = 1, LCK 3 FXY(I) = CONST C RETURN END SUBROUTINE F14 (X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR PROBLEM 12 C DOUBLE PRECISION X, Y(LCK), FXY(LCK), C2PY(1000), C2PX, 1 CONST INTEGER INIT COMMON /INSAV/ C2PY, INIT C CONST = 1.0D0 / DEXP(0.5D0) IF (INIT .GT. 0) GO TO 3 DO 2 I = 1, LCK IF (Y(I) .LE. 0.5D0) GO TO 1 C2PY(I) = - CONST * DEXP (Y(I)) GO TO 2 1 C2PY(I) = 0.0D0 2 CONTINUE INIT = 1 3 CONTINUE C IF (X .LE. 0.5D0) GO TO 4 C2PX = - CONST * DEXP(X) GO TO 5 4 C2PX = 0.0D0 5 CONTINUE DO 6 I = 1, LCK 6 FXY(I) = C2PY(I) + C2PX C RETURN END SUBROUTINE UTRU12 (X, Y, UTRU) C C IN UTRU(1) THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR PROBLEM 10 C IN UTRUE(2) ITS DERIVATIVE WITH RESPECT TO X IS C RETURNED AND IN UTRUE(3) ITS DERIVATIVE WITH RESPECT C TO Y IS RETURNED C DOUBLE PRECISION X, Y, UTRU(3), V, VX, W, WY, CONST C DATA CONST /0.0D0/ IF (CONST. EQ.0.0D0)CONST = 1.0D0/ (2.0D0 * DEXP(0.5D0)) V = X VX = 1.0D0 C IF (Y .LT. 0.5D0) GO TO 1 W = CONST * DEXP(Y) WY = W GO TO 2 1 W = Y WY = 1.0D0 2 CONTINUE C UTRU(1) = V + W UTRU(2) = VX UTRU(3) = WY RETURN END SUBROUTINE UTRU13 (X, Y, UTRU) C C IN UTRU(1) THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR PROBLEM 11 C IN UTRUE(2) ITS DERIVATIVE WITH RESPECT TO X IS C RETURNED AND IN UTRUE(3) ITS DERIVATIVE WITH RESPECT C TO Y IS RETURNED C DOUBLE PRECISION X, Y, UTRU(3), V, VX, W, WY, CONST C DATA CONST /0.0D0/ IF (CONST. EQ.0.0D0)CONST = 1.0D0/ (2.0D0 * DEXP(0.5D0)) W = Y WY = 1.0D0 C IF (X .LT. 0.5D0) GO TO 1 V = CONST * DEXP(X) VX = V GO TO 2 1 V = X VX = 1.0D0 2 CONTINUE C UTRU(1) = V + W UTRU(2) = VX UTRU(3) = WY RETURN END SUBROUTINE UTRU14 (X, Y, UTRU) C C IN UTRU(1) THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR PROBLEM 12 C IN UTRUE(2) ITS DERIVATIVE WITH RESPECT TO X IS C RETURNED AND IN UTRUE(3) ITS DERIVATIVE WITH RESPECT C TO Y IS RETURNED C DOUBLE PRECISION X, Y, UTRU(3), V, VX, W, WY, CONST C DATA CONST /0.0D0/ IF (CONST. EQ.0.0D0)CONST = 1.0D0/ (2.0D0 * DEXP(0.5D0)) C IF (X .LT. 0.5D0) GO TO 1 V = CONST * DEXP(X) VX = V GO TO 2 1 V = X VX = 1.0D0 2 CONTINUE C IF (Y .LT. 0.5D0) GO TO 3 W = CONST * DEXP(Y) WY = W GO TO 4 3 W = Y WY = 1.0D0 4 CONTINUE C UTRU(1) = V + W UTRU(2) = VX UTRU(3) = WY RETURN END DOUBLE PRECISION FUNCTION GAMD6(X) C THIS FUNCTION DEFINES A NEUMANN BOUNDARY FOR PROBLEM 2 DOUBLE PRECISION X GAMD6=3.D0*DEXP(1.D0+X)*X*(1.D0-X) RETURN END DOUBLE PRECISION FUNCTION GAMB7(Y) C THIS FUNCTION DEFINES A DIRICHLET CONDITION FOR PROBLEM 3 DOUBLE PRECISION Y GAMB7=0.75D0*DEXP(0.5D0+Y)*(Y*(1.D0-Y)) RETURN END DOUBLE PRECISION FUNCTION GAMB9(Y) C THIS FUNCTION DEFINES A DIRICHLET BOUNDARY FOR PROBLEMS 5,6,7, AND 9 DOUBLE PRECISION Y GAMB9=DCOS(8.D0*Y*(DATAN(1.0D0))) RETURN END DOUBLE PRECISION FUNCTION GAMB11(Y) C THIS SUBROUTINE DEFINES A NEUMANN CONDITION FOR PROBLEM 9 DOUBLE PRECISION Y GAMB11=-2.0D0*DCOS(8.D0*Y*(DATAN(1.0D0))) RETURN END DOUBLE PRECISION FUNCTION GAMC9(X) C THIS FUNCTION DEFINES A DIRICHLET CONDITION FOR PROBLEM 5 DOUBLE PRECISION X GAMC9=X*X RETURN END DOUBLE PRECISION FUNCTION GAMA12(Y) C THIS FUNCTION DEFINES A DIRICHLET CONDITION FOR PROBLEMS 10 AND 11 DOUBLE PRECISION Y IF (Y .LE. 0.5D0) GO TO 1 GAMA12 = 0.5D0 / DEXP(0.5D0) * DEXP(Y) GO TO 2 1 GAMA12 = Y 2 CONTINUE RETURN END DOUBLE PRECISION FUNCTION GAMB12 (Y) C THIS FUNCTION DEFINES A DIRICHLET CONDITION FOR PROBLEM 10 DOUBLE PRECISION Y, GAMA12 GAMB12 = 1.0D0 + GAMA12(Y) RETURN END DOUBLE PRECISION FUNCTION GAMB14(Y) C THIS FUNCTION DEFINES A DIRICHLET CONDITION FOR PROBLEM 12 DOUBLE PRECISION Y, GAMA12 GAMB14 = GAMA12(Y) + DEXP(0.5D0) / 2.0D0 RETURN END DOUBLE PRECISION FUNCTION GAMC12(X) C THIS SUBROUTINE DEFINES A DIRICHLET CONDITION FOR PROBLEMS 10 AND 11 DOUBLE PRECISION X GAMC12 = X RETURN END DOUBLE PRECISION FUNCTION GAMD12 (X) C THIS SUBROUTINE DEFINES A DIRICHLET CONDITION FOR PROBLEMS 10 AND 11 DOUBLE PRECISION X GAMD12 = X + DEXP(0.5D0) / 2.0D0 RETURN END DOUBLE PRECISION FUNCTION P212(Y) C THIS FUNCTION DEFINES P(Y) FOR PROBLEMS 10 AND 11 DOUBLE PRECISION Y IF (Y .LE. 0.5D0) GO TO 1 P212 = 2.0D0 GO TO 2 1 P212 = 1.0D0 2 CONTINUE RETURN END SUBROUTINE EREST3(K, M, X,IX, N, Y,IY, U, UTRUE,ERROR) INTEGER K, M, N DOUBLE PRECISION X(1), Y(1), U(M,N) EXTERNAL UTRUE C C THIS SUBROUTINE COMPUTES THE MAXIMUM ERROR ON 10,000 C UNIFORMLY SPACED POINTS DEFINED ON THE TENSOR PRODUCT C MESH GIVEN IN X AND Y C C INPUT PARAMETERS C K IS THE ORDER C M IS THE DIMENSION OF BI(X) C X IS THE KNOT SEQUENCE IN THE X DIRECTION C IX IF 0, X IS A NORMAL MESH AND OTHERWISE A PERIODIC MESH C N IS THE DIMENSION OF BJ(X) C Y IS THE KNOT SEQUENCE IN THE Y DIRECTION C IY IF 0, Y IS A NORMAL MESH AND OTHERWISE A PERIODIC ONE C U IS THE M BY N ARRAY OF COEFFICIENTS C UTRUE IS A FUNCTION FOR EVALUTING THE TRUE VALUE AT (X,Y) C AND ITS DERIVATIVES C C OUTPUT PARAMETERS C ERROR ERROR(1) IS THE MAXIMUM ERROR ON 10,000 UNIFOMLY C PLACED POINTS. ERROR(2) IS THE MAXIMUM ERROR IN THE C DERIVATIVE WITH RESPECT TO X AT THESE POINTS AND ERROR(3) C IS THE MAXIMUM ERROR WITH RESPECT TO Y AT THESE POINTS C DOUBLE PRECISION UC(101,101,3), XPTS(101), YPTS(101), * DX, DY, ERROR(3), TERROR,TRUEC(3) C C COMPUTE THE MESH C DX = (X(M+1) - X(K))/100.D0 DY = (Y(N+1) - Y(K))/100.D0 DO 10 I=1,101 XPTS(I) = X(K)+ DBLE(FLOAT (I-1))*DX 10 YPTS(I) = Y(K) + DBLE(FLOAT (I-1))*DY C CALL B2EVAL(K, M, X,IX, N, Y,IY, U, 101, XPTS, 101, YPTS, 1 1, 101,UC,UC(1,1,2),UC(1,1,3)) C C COMPUTE THE MAXIMUM ERROR C ERROR(1)=0.0D0 ERROR(2)=0.0D0 ERROR(3)=0.0D0 DO 20 I=1,101 DO 20 J=1,101 CALL UTRUE(XPTS(I),YPTS(J),TRUEC) DO 18 L=1,3 TERROR=DABS(UC(I,J,L)-TRUEC(L)) IF (TERROR.GT.ERROR(L))ERROR(L)=TERROR 18 CONTINUE 20 CONTINUE RETURN END C-----TOMS TEST.DAT(INPUT TO TEST) 2 2 1 0.D0 1.D0 0.D0 1.D0 2 3 1 0.D0 1.D0 0.D0 1.D0 2 4 1 0.D0 1.D0 0.D0 1.D0 2 5 1 0.D0 1.D0 0.D0 1.D0 4 3 1 0.D0 1.D0 0.D0 1.D0 4 4 1 0.D0 1.D0 0.D0 1.D0 4 5 1 0.D0 1.D0 0.D0 1.D0 6 3 1 0.D0 1.D0 0.D0 1.D0 6 4 1 0.D0 1.D0 0.D0 1.D0 2 2 2 0.D0 1.D0 0.D0 1.D0 2 3 2 0.D0 1.D0 0.D0 1.D0 2 4 2 0.D0 1.D0 0.D0 1.D0 2 5 2 0.D0 1.D0 0.D0 1.D0 4 3 2 0.D0 1.D0 0.D0 1.D0 4 4 2 0.D0 1.D0 0.D0 1.D0 4 5 2 0.D0 1.D0 0.D0 1.D0 6 3 2 0.D0 1.D0 0.D0 1.D0 6 4 2 0.D0 1.D0 0.D0 1.D0 2 2 3 0.D0 0.5D0 0.D0 0.5D0 2 3 3 0.D0 0.5D0 0.D0 0.5D0 2 4 3 0.D0 0.5D0 0.D0 0.5D0 2 5 3 0.D0 0.5D0 0.D0 0.5D0 4 3 3 0.D0 0.5D0 0.D0 0.5D0 4 4 3 0.D0 0.5D0 0.D0 0.5D0 4 5 3 0.D0 0.5D0 0.D0 0.5D0 6 3 3 0.D0 0.5D0 0.D0 0.5D0 6 4 3 0.D0 0.50 0.D0 0.50 2 2 4 0.D0 0.50 0.D0 0.50 2 3 4 0.D0 0.50 0.D0 0.50 2 4 4 0.D0 0.50 0.D0 0.50 2 5 4 0.D0 0.50 0.D0 0.50 4 3 4 0.D0 0.50 0.D0 0.50 4 4 4 0.D0 0.50 0.D0 0.50 4 5 4 0.D0 0.50 0.D0 0.50 6 3 4 0.D0 0.50 0.D0 0.50 6 4 4 0.D0 0.50 0.D0 0.50 2 2 5 0.D0 1.D0 0.D0 1.D0 2 3 5 0.D0 1.D0 0.D0 1.D0 2 4 5 0.D0 1.D0 0.D0 1.D0 2 5 5 0.D0 1.D0 0.D0 1.D0 4 3 5 0.D0 1.D0 0.D0 1.D0 4 4 5 0.D0 1.D0 0.D0 1.D0 4 5 5 0.D0 1.D0 0.D0 1.D0 6 3 5 0.D0 1.D0 0.D0 1.D0 6 4 5 0.D0 1.D0 0.D0 1.D0 2 2 6 0.D0 1.D0 0.D0 1.D0 2 3 6 0.D0 1.D0 0.D0 1.D0 2 4 6 0.D0 1.D0 0.D0 1.D0 2 5 6 0.D0 1.D0 0.D0 1.D0 4 3 6 0.D0 1.D0 0.D0 1.D0 4 4 6 0.D0 1.D0 0.D0 1.D0 4 5 6 0.D0 1.D0 0.D0 1.D0 6 3 6 0.D0 1.D0 0.D0 1.D0 6 4 6 0.D0 1.D0 0.D0 1.D0 2 2 7 0.D0 1.D0 0.D0 1.D0 2 3 7 0.D0 1.D0 0.D0 1.D0 2 4 7 0.D0 1.D0 0.D0 1.D0 2 5 7 0.D0 1.D0 0.D0 1.D0 4 3 7 0.D0 1.D0 0.D0 1.D0 4 4 7 0.D0 1.D0 0.D0 1.D0 4 5 7 0.D0 1.D0 0.D0 1.D0 6 3 7 0.D0 1.D0 0.D0 1.D0 6 4 7 0.D0 1.D0 0.D0 1.D0 2 2 8 0.D0 1.D0 0.D0 1.D0 2 3 8 0.D0 1.D0 0.D0 1.D0 2 4 8 0.D0 1.D0 0.D0 1.D0 2 5 8 0.D0 1.D0 0.D0 1.D0 4 3 8 0.D0 1.D0 0.D0 1.D0 4 4 8 0.D0 1.D0 0.D0 1.D0 4 5 8 0.D0 1.D0 0.D0 1.D0 6 3 8 0.D0 1.D0 0.D0 1.D0 6 4 8 0.D0 1.D0 0.D0 1.D0 2 2 9 0.D0 1.D0 0.D0 1.D0 2 3 9 0.D0 1.D0 0.D0 1.D0 2 4 9 0.D0 1.D0 0.D0 1.D0 2 5 9 0.D0 1.D0 0.D0 1.D0 4 3 9 0.D0 1.D0 0.D0 1.D0 4 4 9 0.D0 1.D0 0.D0 1.D0 4 5 9 0.D0 1.D0 0.D0 1.D0 6 3 9 0.D0 1.D0 0.D0 1.D0 6 4 9 0.D0 1.D0 0.D0 1.D0 2 2 10 0.D0 1.D0 0.D0 1.D0 2 3 10 0.D0 1.D0 0.D0 1.D0 2 4 10 0.D0 1.D0 0.D0 1.D0 2 5 10 0.D0 1.D0 0.D0 1.D0 4 3 10 0.D0 1.D0 0.D0 1.D0 4 4 10 0.D0 1.D0 0.D0 1.D0 4 5 10 0.D0 1.D0 0.D0 1.D0 6 3 10 0.D0 1.D0 0.D0 1.D0 6 4 10 0.D0 1.D0 0.D0 1.D0 2 2 11 0.D0 1.D0 0.D0 1.D0 2 3 11 0.D0 1.D0 0.D0 1.D0 2 4 11 0.D0 1.D0 0.D0 1.D0 2 5 11 0.D0 1.D0 0.D0 1.D0 4 3 11 0.D0 1.D0 0.D0 1.D0 4 4 11 0.D0 1.D0 0.D0 1.D0 4 5 11 0.D0 1.D0 0.D0 1.D0 6 3 11 0.D0 1.D0 0.D0 1.D0 6 4 11 0.D0 1.D0 0.D0 1.D0 2 2 12 0.D0 1.D0 0.D0 1.D0 2 3 12 0.D0 1.D0 0.D0 1.D0 2 4 12 0.D0 1.D0 0.D0 1.D0 2 5 12 0.D0 1.D0 0.D0 1.D0 4 3 12 0.D0 1.D0 0.D0 1.D0 4 4 12 0.D0 1.D0 0.D0 1.D0 4 5 12 0.D0 1.D0 0.D0 1.D0 6 3 12 0.D0 1.D0 0.D0 1.D0 6 4 12 0.D0 1.D0 0.D0 1.D0 0 0 0 C-----TOMS TEST.OUT (OUTPUT FROM TEST) This is the output file produced on a Vax 750. The problem number corresponds to the test problems( not the example problems). The order is the order of the approximation. If the mesh size is given as p then 2**p points have been put on a side. The error in the solution and the derivatives is the maximum determined at 10000 uniformly spaced points. The error in the solution should decrease by a factor of 4 when the order is 2, by a factor of 16 when the order is 4, and a factor of 64 when the order is 6. The error in the derivatives shoiuld decrease be a factor of 2 when the order is 2, a factor of 8 when the order is 4, and a factor of 32 when the order is 16. PROBLEM ORDER MESH SIZE ERROR ERROR ERROR SOLUTION DERIVATIVE DERIVATIVE IN X IN Y 1 2 5 .738E-01 .140E+01 .140E+01 1 2 9 .222E-01 .774E+00 .774E+00 1 2 17 .619E-02 .415E+00 .415E+00 1 2 33 .150E-02 .215E+00 .215E+00 1 4 9 .254E-04 .109E-02 .109E-02 1 4 17 .150E-05 .133E-03 .133E-03 1 4 33 .923E-07 .163E-04 .163E-04 1 6 9 .260E-07 .729E-06 .729E-06 1 6 17 .347E-09 .254E-07 .254E-07 2 2 5 .201E+00 .143E+01 .143E+01 2 2 9 .497E-01 .771E+00 .771E+00 2 2 17 .124E-01 .415E+00 .415E+00 2 2 33 .309E-02 .216E+00 .216E+00 2 4 9 .254E-04 .109E-02 .109E-02 2 4 17 .150E-05 .133E-03 .133E-03 2 4 33 .923E-07 .163E-04 .163E-04 2 6 9 .260E-07 .729E-06 .729E-06 2 6 17 .347E-09 .254E-07 .254E-07 3 2 5 .712E-02 .204E+00 .204E+00 3 2 9 .206E-02 .107E+00 .107E+00 3 2 17 .548E-03 .546E-01 .546E-01 3 2 33 .125E-03 .276E-01 .276E-01 3 4 9 .792E-06 .648E-04 .648E-04 3 4 17 .563E-07 .784E-05 .784E-05 3 4 33 .342E-08 .965E-06 .965E-06 3 6 9 .232E-09 .370E-07 .370E-07 3 6 17 .394E-11 .120E-08 .120E-08 4 2 5 .164E-01 .188E+00 .188E+00 4 2 9 .473E-02 .987E-01 .987E-01 4 2 17 .127E-02 .525E-01 .525E-01 4 2 33 .312E-03 .270E-01 .270E-01 4 4 9 .796E-06 .542E-04 .542E-04 4 4 17 .565E-07 .678E-05 .678E-05 4 4 33 .339E-08 .847E-06 .847E-06 4 6 9 .232E-09 .117E-07 .117E-07 4 6 17 .394E-11 .390E-09 .390E-09 5 2 5 .222E+00 .155E+01 .489E+01 5 2 9 .540E-01 .464E+00 .247E+01 5 2 17 .149E-01 .183E+00 .123E+01 5 2 33 .357E-02 .867E-01 .617E+00 5 4 9 .572E-03 .283E-02 .397E-01 5 4 17 .341E-04 .200E-03 .406E-02 5 4 33 .208E-05 .208E-04 .483E-03 5 6 9 .997E-05 .302E-03 .117E-02 5 6 17 .131E-06 .837E-05 .360E-04 6 2 5 .222E+00 .140E+01 .489E+01 6 2 9 .527E-01 .458E+00 .247E+01 6 2 17 .129E-01 .156E+00 .123E+01 6 2 33 .322E-02 .556E-01 .617E+00 6 4 9 .610E-03 .127E-02 .260E-01 6 4 17 .343E-04 .692E-04 .309E-02 6 4 33 .208E-05 .418E-05 .378E-03 6 6 9 .105E-04 .211E-04 .469E-03 6 6 17 .131E-06 .263E-06 .126E-04 7 2 5 .222E+00 .489E+01 .149E+01 7 2 9 .527E-01 .247E+01 .461E+00 7 2 17 .129E-01 .123E+01 .155E+00 7 2 33 .322E-02 .617E+00 .552E-01 7 4 9 .610E-03 .260E-01 .127E-02 7 4 17 .343E-04 .309E-02 .692E-04 7 4 33 .208E-05 .378E-03 .418E-05 7 6 9 .105E-04 .469E-03 .211E-04 7 6 17 .131E-06 .126E-04 .263E-06 8 2 5 .260E+00 .579E+00 .366E+01 8 2 9 .830E-01 .230E+00 .231E+01 8 2 17 .220E-01 .892E-01 .121E+01 8 2 33 .559E-02 .379E-01 .614E+00 8 4 9 .603E-03 .122E-02 .259E-01 8 4 17 .342E-04 .686E-04 .309E-02 8 4 33 .208E-05 .417E-05 .378E-03 8 6 9 .105E-04 .210E-04 .469E-03 8 6 17 .131E-06 .263E-06 .126E-04 9 2 5 .211E+00 .400E+01 .585E+00 9 2 9 .678E-01 .235E+01 .221E+00 9 2 17 .180E-01 .122E+01 .887E-01 9 2 33 .455E-02 .615E+00 .374E-01 9 4 9 .603E-03 .259E-01 .122E-02 9 4 17 .342E-04 .309E-02 .686E-04 9 4 33 .208E-05 .378E-03 .417E-05 9 6 9 .105E-04 .469E-03 .210E-04 9 6 17 .131E-06 .126E-04 .263E-06 10 2 5 .536E-02 .122E-01 .945E-01 10 2 9 .147E-02 .527E-02 .487E-01 10 2 17 .383E-03 .264E-02 .250E-01 10 2 33 .900E-04 .133E-02 .127E-01 10 4 9 .262E-06 .116E-05 .184E-04 10 4 17 .161E-07 .997E-07 .208E-05 10 4 33 .972E-09 .113E-07 .255E-06 10 6 9 .106E-09 .400E-08 .130E-07 10 6 17 .177E-11 .128E-09 .445E-09 11 2 5 .536E-02 .945E-01 .122E-01 11 2 9 .147E-02 .487E-01 .527E-02 11 2 17 .383E-03 .250E-01 .264E-02 11 2 33 .900E-04 .127E-01 .133E-02 11 4 9 .262E-06 .184E-04 .116E-05 11 4 17 .161E-07 .208E-05 .997E-07 11 4 33 .972E-09 .255E-06 .113E-07 11 6 9 .106E-09 .130E-07 .400E-08 11 6 17 .177E-11 .445E-09 .128E-09 12 2 5 .793E-02 .100E+00 .100E+00 12 2 9 .207E-02 .510E-01 .510E-01 12 2 17 .526E-03 .260E-01 .260E-01 12 2 33 .125E-03 .132E-01 .132E-01 12 4 9 .497E-06 .195E-04 .195E-04 12 4 17 .322E-07 .218E-05 .218E-05 12 4 33 .190E-08 .266E-06 .266E-06 12 6 9 .211E-09 .150E-07 .150E-07 12 6 17 .353E-11 .512E-09 .512E-09 C-----TOMS EXAMPLE (MAIN PROGRAM) C SERRG2D EXAMPLES C C THIS PROGRAM SOLVES EXAMPLES 1 THRU 5 IN THE PAPER C "A PROGRAM FOR SOLVING SEPARABLE ELLIPTIC EQUATIONS" C BY LINDA KAUFMAN AND DANIEL D. WARNER C DOUBLE PRECISION X(129), Y(129), V(4900), ERROR, ZERO, ONE, 1 FOUR, A, B, C, D, VALKNT, ERR1, ERR2, 1 F1, UTRUE1, 1 F2, UTRUE2, 1 F3, UTRUE3, GAMB3, 1 F4, UTRUE4, P14, R4, GAMB4, 1 F5, UTRUE5, P25, GAMA5, GAMB5, DX EXTERNAL ZERO, ONE, FOUR, 1 F1, UTRUE1, 1 F2, UTRUE2, 1 F3, UTRUE3, GAMB3, 1 F4, UTRUE4, P14, R4, GAMB4, 1 F5, UTRUE5, P25, GAMA5, GAMB5 C INTEGER INIT, INUNIT, OUTUNT DOUBLE PRECISION SAVE(1000) COMMON /INSAV/ SAVE, INIT C COMMON /CSTAK/ DSTAK(100000) DOUBLE PRECISION DSTAK C CALL ISTKIN(100000, 4) OUTUNT = I1MACH(2) INUNIT = I1MACH(1) WRITE(OUTUNT,5) 5 FORMAT(8H EXAMPLE,4X,5HORDER,4X,4HMESH,3X, 1 13HUNIFORM ERROR,4X,12HGRADED ERROR) C 10 READ (INUNIT, 20) K, KK, IEXMPL 20 FORMAT (3I5) IF (K .EQ. 0) STOP MESH = 2**(KK) + 1 INIT = 0 GO TO (50,60,70,80,90), IEXMPL C C EXAMPLE 1 POISSON WITH DIRICHLET BOUNDARY UNIFORM MESH C 50 A = 0.0D0 B = 1.0D0 C = 0.0D0 D = 1.0D0 C CALL UMESH (A, B, MESH, K, X, M) CALL UMESH (C, D, MESH, K, Y, N) C C SOLVE THE PDE CALL SERRG2(K, M, X, N, Y, 1 ONE, ZERO, ZERO, ONE, ZERO, F1, V, 1 1, 0.0D0, ZERO, 1, 0.0D0, ZERO, 1 1, 0.0D0, ZERO, 1, 0.0D0, ZERO) C C DETERMINE THE MAXIMUM ERROR AT 10000 UNIFORMLY SPACED C POINTS C CALL ERREST(K, M, X,0, N, Y,0, V, UTRUE1, ERROR) C WRITE (OUTUNT,57) IEXMPL, K, MESH, ERROR 57 FORMAT(1H ,I4,I11,I8,D15.3,D17.3) GOTO 10 C C EXAMPLE 2 POISSON, DIRICHLET,SINGULAR SECOND DERIVATIVES C 60 A = 0.0D0 B = 1.0D0 C = 0.0D0 D = 1.0D0 N = MESH + K - 2 M = N C C GENERATE GRADED MESH C DO 64 I = 1, K X(I) = A Y(I) = C NPI=N+I X(NPI) = B Y(NPI) = D 64 CONTINUE DX=1.D0/(DBLE(FLOAT(MESH))-1.D0) DO 66 I = 1, MESH IK1 = I + K - 1 X(IK1) = (DBLE(FLOAT(I-1))*DX)**3 Y(IK1) = X(IK1) 66 CONTINUE C C SOLVE THE PDE CALL SERRG2(K, M, X, N, Y, 1 ONE, ZERO, ZERO, ONE, ZERO, F2, V, 1 1, 0.0D0, ZERO, 1, 0.0D0, ZERO, 1 1, 0.0D0, ZERO, 1, 0.0D0, ZERO) C C DETERMINE THE MAXIMUM ERROR AT 10000 UNIFORMLY SPACED C POINTS C CALL ERREST(K, M, X,0, N, Y,0, V, UTRUE2, ERROR) ERR2 = ERROR C C GENERATE UNIFORM MESH C INIT = 0 CALL UMESH (A, B, MESH, K, X, M) CALL UMESH (C, D, MESH, K, Y, N) C C SOLVE THE PDE CALL SERRG2(K, M, X, N, Y, 1 ONE, ZERO, ZERO, ONE, ZERO, F2, V, 1 1, 0.0D0, ZERO, 1, 0.0D0, ZERO, 1 1, 0.0D0, ZERO, 1, 0.0D0, ZERO) C C DETERMINE THE MAXIMUM ERROR AT 10000 UNIFORMLY SPACED C POINTS C CALL ERREST(K, M, X,0, N, Y,0, V, UTRUE2, ERROR) ERR1 = ERROR C WRITE (OUTUNT,57) IEXMPL, K, MESH, ERR1, ERR2 GO TO 10 C C EXAMPLE 3 CONSTANT COEFFICIENT. DIRICHELET, NEUMANN IN X PERIODIC IN Y C UNIFORM MESH C 70 A = 0.0D0 B = 1.0D0 C = 0.0D0 D = 1.0D0 C CALL UMESH (A, B, MESH, K, X, M) DO 71 I=1, MESH 71 Y(I)=DBLE(FLOAT(I-1))/DBLE(FLOAT(MESH-1)) N = MESH - 1 C C SOLVE THE PROBLEM ON A UNIFORM MESH CALL SERRG2(K, M, X, N, Y, 1 ONE, ZERO, ZERO, ONE, FOUR, F3, V, 1 1, 0.0D0, ZERO, 3, 0.0D0, GAMB3, 1 4, 0.0D0, ZERO, 1, 0.0D0, ZERO) C C DETERMINE THE MAXIMUM ERROR AT 10000 UNIFORMLY SPACED C POINTS C CALL ERREST(K, M, X,0, N, Y,1, V, UTRUE3, ERROR) C WRITE (OUTUNT,57) IEXMPL, K, MESH, ERROR GO TO 10 C C EXAMPLE 4 POISSON IN POLAR COOR. DIRICHLET AND NEUMANN C 80 A = 0.0D0 B = 1.0D0 C = 0.0D0 D = 2.0D0 * DATAN(1.0D0) C CALL UMESH (A, B, MESH, K, X, M) CALL UMESH (C, D, MESH, K, Y, N) C C SOLVE THE PDE ON UNIFORM MESHES CALL SERRG2(K, M, X, N, Y, 1 P14, R4, ZERO, ONE, ZERO, F4, V, 1 1, 0.0D0, ZERO, 1, 0.0D0, GAMB4, 1 3, 0.0D0, ZERO, 3, 0.0D0, ZERO ) C C DETERMINE THE MAXIMUM ERROR AT 10000 UNIFORMLY SPACED C POINTS C CALL ERREST(K, M, X, 0, N,Y,0,V, UTRUE4,ERROR) C WRITE (OUTUNT,57) IEXMPL, K, MESH, ERROR GOTO 10 C C EXAMPLE 5 ELECTROSTATIC POTENTIAL ON A WASHER C 90 A = 1.0D0 B = 2.0D0 C = 0.0D0 D = 4.0D0 * DATAN(1.0D0) C CALL UMESH (A, B, MESH, K, X, M) C C MAKE A MESH WITH MULTIPLICITY K-1 AT PI/2 BUT C UNIFORM IN BOTH SEGMENTS C VALKNT = DATAN(1.0D0) CALL BKMSH(C, D, MESH, K, VALKNT, Y, N) C SOLVE THE PDE CALL SERRG2(K, M, X, N, Y, 1 P14, R4, ZERO, P25, ZERO, F5, V, 1 1, 0.0D0, GAMA5, 1, 0.0D0, GAMB5, 1 3, 0.0D0, ZERO, 3, 0.0D0, ZERO) C C DETERMINE THE MAXIMUM ERROR AT 10000 UNIFORMLY SPACED C POINTS C CALL ERREST(K, M, X, 0, N,Y,0,V, UTRUE5,ERROR) C WRITE (OUTUNT,57) IEXMPL, K, MESH, ERROR GOTO 10 C END SUBROUTINE ERREST(K, M, X,IX, N, Y,IY, U, UTRUE,ERROR) INTEGER K, M, N DOUBLE PRECISION X(1), Y(1), U(M,N) EXTERNAL UTRUE C C THIS SUBROUTINE COMPUTES THE MAXIMUM ERROR ON 10,000 C UNIFORMLY SPACED POINTS DEFINED ON THE TENSOR PRODUCT C MESH GIVEN IN X AND Y C C INPUT PARAMETERS C K IS THE ORDER C M IS THE DIMENSION OF BI(X) C X IS THE KNOT SEQUENCE IN THE X DIRECTION C IX IF 0, X IS A NORMAL MESH AND OTHERWISE A PERIODIC MESH C N IS THE DIMENSION OF BJ(X) C Y IS THE KNOT SEQUENCE IN THE Y DIRECTION C IY IF 0, Y IS A NORMAL MESH AND OTHERWISE A PERIODIC ONE C U IS THE M BY N ARRAY OF COEFFICIENTS C UTRUE IS A FUNCTION FOR EVALUTING THE TRUE VALUE AT (X,Y) C C OUTPUT PARAMETERS C ERROR THE MAXIMUM ERROR ON 10,000 UNIFOMLY PLACED POINTS C DOUBLE PRECISION UC(101,101), XPTS(101), YPTS(101), 1 DX, DY, ERROR, TERROR,TRUEC,UD1,UD2 C C COMPUTE THE MESH C DX = (X(M+1) - X(K))/100.D0 DY = (Y(N+1) - Y(K))/100.D0 DO 10 I=1,101 XPTS(I) = X(K)+ DBLE(FLOAT (I-1))*DX 10 YPTS(I) = Y(K) + DBLE(FLOAT (I-1))*DY C CALL B2EVAL(K, M, X,IX, N, Y,IY, U, 101, XPTS, 101, YPTS, 1 0, 101,UC,UD1,UD2) C C COMPUTE THE MAXIMUM ERROR C ERROR=0.0D0 DO 20 I=1,101 DO 20 J=1,101 CALL UTRUE(XPTS(I),YPTS(J),TRUEC) TERROR=DABS(UC(I,J)-TRUEC) IF (TERROR.GT.ERROR)ERROR=TERROR 20 CONTINUE RETURN END SUBROUTINE UMESH(A, B, NX, K, T, N) C THIS SUBROUTINE COMPUTES A UNIFORM MESH WITH NX-1 INTERVALS C WITH THE ENDPOINTS AT A AND B HAVING MULTIPLICITY K C C INPUT PARAMETERS C C A THE VALUE OF THE FIRST POINT IN THE MESH C B THE VALUE OF THE LAST POINT OF THE MESH C NX THE NUMBER OF INTERVALS -1 IN THE UNIFORM MESH C K NUMBER OF MULTIPLE KNOTS AT THE END POINTS C C OUTPUT PARAMETERS C C T THE KNOT SEQUENCE C N K+NX-2, THE NUMBER OF B-SPLINES INTEGER NX, K, N DOUBLE PRECISION A, B, T(1), DX C C LEFT HAND ENDPOINTS DO 10 I=1,K 10 T(I) = A DX = (B - A)/DBLE(FLOAT(NX-1)) C COMPUTE THE INTERIOR KNOTS DO 20 I=1,NX IK1 = I + K - 1 20 T(IK1) = DBLE(FLOAT(I-1))*DX + A N = K + NX - 2 C RIGHT HAND ENDPOINTS DO 30 I=1,K NPI = N+I 30 T(NPI) = B RETURN END SUBROUTINE BKMSH(A, B, NX, K, BRK, T, N) C THIS SUBROUTINE CREATES A MESH WHICH HAS ENDPOINTS C A AND B AND A KNOT OF MULTIPLICITY K-1 AT POINT BRK C AND HAS UNIFORM INTERVALS BEFORE AND AFTER BRK WITH C NX KNOTS SPREAD PROPROTIONATELY BEFORE AND AFTER BRK C C INPUT PARAMETERS C C A THE VALUE OF THE FIRST POINT IN THE MESH C B THE VALUE OF THE LAST POINT OF THE MESH C NX THE NUMBER OF POINTS IN THE UNIFORM MESH C K NUMBER OF MULTIPLE KNOTS AT THE END POINTS C BRK THE POINT OF THE MULTIPLE KNOT C OUTPUT PARAMETERS C T THE KNOT SEQUENCE C N T HAS N+K ENTRIES INTEGER NX, K, N DOUBLE PRECISION A, B, T(1), DX DOUBLE PRECISION BRK, DX2 C C LEFT HAND ENDPOINTS DO 10 I=1,K 10 T(I) = A C NX1=(BRK-A)/(B-A)*DBLE(FLOAT(NX))+1.D0 NX2 = NX-NX1+1 IF (NX1.LT.1) GO TO 30 C C DO INTERIOR POINTS IN THE FIRST INTERVAL C DX=(BRK-A)/DBLE(FLOAT(NX1-1)) DO 20 I=1,NX1 IK=I+K-1 T(IK)=DBLE(FLOAT(I-1))*DX+A 20 CONTINUE C C DO THE BREAKPOINT C 30 N=NX1+K-2 KM1= K-1 DO 40 I=1,KM1 NPI=N+I 40 T(NPI)=BRK C C DO THE POINTS IN THE SECOND INTERVAL C IF (NX2.LT.1) GO TO 60 DX2= (B-BRK)/DBLE(FLOAT(NX2-1)) NMK2=N+K-2 DO 50 I=1,NX2 IK2 = NMK2+I T(IK2)=BRK+DBLE(FLOAT (I-1))*DX2 50 CONTINUE C C DO THE RIGHT HAND ENDPOINT 60 N=NMK2+NX2-1 DO 70 I=1,K NPI=N+I T(NPI)=B 70 CONTINUE RETURN END DOUBLE PRECISION FUNCTION ZERO(X) C C THIS FUNCTION RETURNS THE VALUE 0 REGARDLESS C OF THE VALUE OF THE PARAMETER X C DOUBLE PRECISION X ZERO = 0.0D0 RETURN END DOUBLE PRECISION FUNCTION ONE(X) C C THIS FUNCTION RETURNS THE VALUE 1 REGARDLESS C OF THE VALUE OF THE PARAMETER X C DOUBLE PRECISION X ONE = 1.0D0 RETURN END DOUBLE PRECISION FUNCTION FOUR(X) C C THIS FUNCTION RETURNS THE VALUE 4 REGARDLESS C OF THE VALUE OF THE PARAMETER X C DOUBLE PRECISION X FOUR=4.0D0 RETURN END SUBROUTINE F1(X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR EXAMPLE 1 C DOUBLE PRECISION X, Y(LCK), FXY(LCK), T1, T2, T3 DOUBLE PRECISION SAVE(1000) INTEGER INIT COMMON /INSAV/ SAVE, INIT C IF (INIT.GT.0)GO TO 3 DO 2 I=1,LCK 2 SAVE(I)=DEXP(Y(I))*Y(I) INIT=1 3 CONTINUE C T1 = -3.0D0 * DEXP(X)*X T2 = X + 3.0D0 T3 = X - 1.0D0 DO 1 I=1,LCK 1 FXY(I) = T1 * SAVE(I) * ((Y(I) - 1.0D0) * T2 + 1 T3 * (Y(I) + 3.0D0)) RETURN END SUBROUTINE UTRUE1(X, Y,UTRU) C C IN UTRU THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR EXAMPLE 1 C DOUBLE PRECISION X, Y,UTRU C UTRU= 3.0D0 * DEXP(X + Y) * X * (X - 1.0D0) * Y * (Y - 1.0D0) RETURN END SUBROUTINE F2(X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR EXAMPLE 2 C DOUBLE PRECISION X, Y(LCK), FXY(LCK), T1, T2, SQRTX, SQRTY DOUBLE PRECISION SAVE(1000) INTEGER INIT COMMON /INSAV/ SAVE, INIT C IF (INIT.GT.0)GO TO 3 DO 2 I=1,LCK SQRTY = DSQRT(Y(I)) SAVE(I) = Y(I) * (SQRTY - 1D0) LCKPI=LCK+I 2 SAVE(LCKPI) = (.75D0) / SQRTY INIT=1 3 CONTINUE C SQRTX = DSQRT(X) T1 = .75D0 / SQRTX T2 = X * (SQRTX - 1D0) DO 1 I=1,LCK LCKPI= LCK+I 1 FXY(I) = - (T1 * SAVE(I) + T2 * SAVE(LCKPI)) RETURN END SUBROUTINE UTRUE2(X, Y,UTRU) C C IN UTRU THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR EXAMPLE 2 C DOUBLE PRECISION X, Y,UTRU, TX0, TY0, SQRTX, SQRTY C SQRTX = DSQRT(X) SQRTY = DSQRT(Y) TX0 = X * (SQRTX - 1D0) TY0 = Y * (SQRTY - 1D0) UTRU= TX0 * TY0 RETURN END SUBROUTINE F3(X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR EXAMPLE 3 C DOUBLE PRECISION X, Y(LCK), FXY(LCK) DOUBLE PRECISION TWOPI, PI2, X2, C2PY(1000),PI2P1,PIO4 INTEGER INIT COMMON /INSAV/ C2PY, INIT DATA PIO4 /0.0D0/ C IF (PIO4.EQ.0.0D0) PIO4=DATAN(1.0D0) PI2 = 16.D0*PIO4**2 PI2P1=PI2+1.D0 IF (INIT.GT.0) GO TO 3 TWOPI=8.D0*PIO4 DO 2 I=1,LCK 2 C2PY(I)=DCOS(TWOPI*Y(I)) INIT=1 3 CONTINUE C X2 = 2.D0-4.D0*X*X*PI2P1 DO 1 I=1,LCK 1 FXY(I) = -C2PY(I)*X2 RETURN END SUBROUTINE UTRUE3(X, Y,UTRU) C C IN UTRU THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR EXAMPLE 3 C DOUBLE PRECISION X, Y,UTRU DOUBLE PRECISION TWOPI DATA TWOPI/0.0D0/ C IF (TWOPI .EQ. 0.0D0) TWOPI = 8.D0*DATAN(1.0D0) UTRU = X*X*(DCOS(TWOPI*Y)) RETURN END DOUBLE PRECISION FUNCTION GAMB3(Y) C C THIS SUBROUTINE RETURNS THE NEUMANN CONDITION FOR C EXAMPLE 3 C DOUBLE PRECISION Y DOUBLE PRECISION TWOPI DATA TWOPI/0.0D0/ IF (TWOPI .EQ. 0.0D0) TWOPI = 8.D0*DATAN(1.0D0) GAMB3=-2.0D0*DCOS(Y*TWOPI) RETURN END SUBROUTINE F4(X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR EXAMPLE 4 C DOUBLE PRECISION X, Y(LCK), FXY(LCK), FF C FF = - 16.D0 * X ** 4 DO 1 I=1,LCK 1 FXY(I) = FF RETURN END SUBROUTINE UTRUE4(X, Y,UTRU) C C IN UTRU THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR EXAMPLE 4 C DOUBLE PRECISION X, Y,UTRU, COS4Y C COS4Y = DCOS(4D0*Y) UTRU= (1D0 - COS4Y) * X ** 4 RETURN END DOUBLE PRECISION FUNCTION P14 (X) C THIS SUBROUTINE DEFINES THE FUNCTION P1(X) FOR EXAMPLE 4 DOUBLE PRECISION X P14 = X * X RETURN END DOUBLE PRECISION FUNCTION R4 (X) C THIS SUBROUTINE DEFINES THE FUNCTION R(X) FOR EXAMPLE 4 DOUBLE PRECISION X R4 = X RETURN END DOUBLE PRECISION FUNCTION GAMB4(X) C THIS SUBROUTINE DEFINES A BOUNDARY FUNCTION FOR EXAMPLE 4 DOUBLE PRECISION X GAMB4 = 1D0 - DCOS(4D0 * X) RETURN END SUBROUTINE F5 (X, Y, FXY, LCK) C C THIS SUBROUTINE RETURNS THE RIGHT HAND C SIDE F(X,Y(I)), I = 1,....LCK IN FXY(I), FOR I=1,2,...LCK C FOR EXAMPLE 5 C DOUBLE PRECISION X, Y(LCK), FXY(LCK), SQRT2, C1, C2, PIO4, 1 SAVE(1000), COSY, COS2Y INTEGER INIT COMMON /INSAV/ SAVE, INIT DATA PIO4 /0.0D0/ DATA SQRT2 /0.0D0/ C IF (SQRT2.EQ.0.0D0)SQRT2=DSQRT(2D0) IF (PIO4.EQ.0.0D0)PIO4=DATAN(1.D0) C1 = (4D0 + SQRT2) / 2D0 C2 = SQRT2 / 2D0 C IF (INIT .GT. 0) GO TO 2 DO 1 I = 1, LCK COSY = DCOS(Y(I)) COS2Y = DCOS(2D0*Y(I)) SAVE(I) = (C1 - 1D0) * COSY IPLCK=I+LCK SAVE(IPLCK) = (4D0 * C2 - 1D0) * COS2Y + (C2 - 1D0) * COSY 1 CONTINUE INIT = 1 2 CONTINUE C DO 4 I = 1, LCK IF (Y(I) .LE. PIO4) GO TO 3 IPLCK=I+LCK FXY(I) = SAVE(IPLCK) * X GO TO 4 3 FXY(I) = SAVE(I) * X 4 CONTINUE C RETURN END SUBROUTINE UTRUE5 (X, Y, UTRU) C C IN UTRU THIS SUBROUTINE RETURNS THE TRUE C SOLUTION AT (X,Y) FOR EXAMPLE 5 C DOUBLE PRECISION X, Y, UTRU, PIO4, COSY, COS2Y, * Y2 DATA PIO4 /0.0D0/ C IF (PIO4.EQ.0.0D0)PIO4=DATAN(1.D0) COSY = DCOS(Y) C IF (Y .GT. PIO4) GO TO 1 UTRU = X * COSY GO TO 2 1 Y2 = 2D0 * Y COS2Y = DCOS(Y2) UTRU = X * (COSY + COS2Y) 2 RETURN END DOUBLE PRECISION FUNCTION P25 (X) C THIS SUBROUTINE DEFINES P(X) FOR EXAMPLE 5 DOUBLE PRECISION X, SQRT2, PIO4 DATA PIO4 /0.0D0/ DATA SQRT2 /0.0D0/ C IF (PIO4.EQ.0.0D0)PIO4=DATAN(1.D0) IF (SQRT2.EQ.0.0D0)SQRT2=DSQRT(2D0) IF (X .GT. PIO4) GO TO 1 P25 = (4D0 + SQRT2) / 2D0 GO TO 2 1 P25 = SQRT2 / 2D0 2 RETURN END DOUBLE PRECISION FUNCTION GAMA5 (X) C THIS SUBROUTINE DEFINES A BOUNDARY CONDITION FOR EXAMPLE 5 DOUBLE PRECISION X, PIO4 DATA PIO4 /0.0D0/ C IF (PIO4.EQ.0.0D0)PIO4=DATAN(1.D0) IF (X .GT. PIO4) GO TO 1 GAMA5 = DCOS(X) GO TO 2 1 GAMA5 = DCOS(X) + DCOS(2D0*X) 2 RETURN END DOUBLE PRECISION FUNCTION GAMB5 (X) C THIS SUBROUTINE DEFINES A BOUNDARY CONDITION FOR EXAMPLE 5 DOUBLE PRECISION X, GAMA5 GAMB5 = 2D0 * GAMA5(X) RETURN END C-----TOMS EXAMPLE.DAT (INPUT TO EXAMPLE) 2 3 1 2 4 1 2 5 1 2 6 1 4 3 1 4 4 1 4 5 1 6 3 1 6 4 1 6 5 1 2 3 2 2 4 2 2 5 2 2 6 2 4 3 2 4 4 2 4 5 2 6 3 2 6 4 2 6 5 2 2 3 3 2 4 3 2 5 3 2 6 3 4 3 3 4 4 3 4 5 3 6 3 3 6 4 3 6 5 3 2 3 4 2 4 4 2 5 4 2 6 4 4 3 4 4 4 4 4 5 4 6 3 4 6 4 4 6 5 4 2 3 5 2 4 5 2 5 5 2 6 5 4 3 5 4 4 5 4 5 5 6 3 5 6 4 5 6 5 5 0 0 0 C-----TOMS EXAMPLE.OUT (OUTPUT FROM EXAMPLE) This is the output file produced on a Vax 750. The problem number corresponds to the examples in the TOMS papers. The order is the order of the approximation. If the mesh size is given as p then 2**p points have been put on a side. The error in the solution for a uniform mesh is the maximum determined at 10000 uniformly spaced points. The error in the solution should decrease by a factor of 4 when the order is 2, by a factor of 16 when the order is 4, and a factor of 64 when the order is 6. For example 2 a graded mesh was also tried. EXAMPLE ORDER MESH UNIFORM ERROR GRADED ERROR 1 2 9 .222E-01 1 2 17 .619E-02 1 2 33 .150E-02 1 2 65 .387E-03 1 4 9 .254E-04 1 4 17 .150E-05 1 4 33 .923E-07 1 6 9 .260E-07 1 6 17 .347E-09 1 6 33 .514E-11 2 2 9 .896E-03 .153E-02 2 2 17 .319E-03 .422E-03 2 2 33 .109E-03 .111E-03 2 2 65 .341E-04 .287E-04 2 4 9 .572E-04 .418E-05 2 4 17 .194E-04 .242E-06 2 4 33 .694E-05 .147E-07 2 6 9 .177E-04 .186E-06 2 6 17 .609E-05 .375E-08 2 6 33 .215E-05 .780E-10 3 2 9 .830E-01 3 2 17 .220E-01 3 2 33 .559E-02 3 2 65 .140E-02 3 4 9 .603E-03 3 4 17 .342E-04 3 4 33 .208E-05 3 6 9 .105E-04 3 6 17 .131E-06 3 6 33 .193E-08 4 2 9 .712E-01 4 2 17 .203E-01 4 2 33 .538E-02 4 2 65 .138E-02 4 4 9 .572E-03 4 4 17 .341E-04 4 4 33 .208E-05 4 6 9 .997E-05 4 6 17 .138E-06 4 6 33 .210E-08 5 2 9 .106E+00 5 2 17 .259E-01 5 2 33 .644E-02 5 2 65 .161E-02 5 4 9 .122E-02 5 4 17 .686E-04 5 4 33 .417E-05 5 6 9 .206E-04 5 6 17 .271E-06 5 6 33 .413E-08