* D03PRF Example Program Text * Mark 16 Release. NAG Copyright 1993. * .. Parameters .. INTEGER NOUT PARAMETER (NOUT=6) INTEGER NPDE, NPTS, NCODE, NXI, NXFIX, NLEFT, NEQN, + LISAVE, NWKRES, LENODE, LRSAVE, INTPTS, ITYPE PARAMETER (NPDE=2,NPTS=61,NCODE=0,NXI=0,NXFIX=0,NLEFT=1, + NEQN=NPDE*NPTS+NCODE,LISAVE=25+NXFIX, + NWKRES=NPDE*(NPTS+21+3*NPDE)+7*NPTS+NXFIX+3, + LENODE=11*NEQN+50,LRSAVE=NEQN*NEQN+NEQN+NWKRES+ + LENODE,INTPTS=5,ITYPE=1) * .. Scalars in Common .. DOUBLE PRECISION P * .. Local Scalars .. DOUBLE PRECISION CON, DXMESH, TOUT, TRMESH, TS, XRATIO, XX INTEGER I, IFAIL, IND, IPMINF, IT, ITASK, ITOL, ITRACE, + NRMESH LOGICAL REMESH, THETA CHARACTER LAOPT, NORM * .. Local Arrays .. DOUBLE PRECISION ALGOPT(30), ATOL(1), RSAVE(LRSAVE), RTOL(1), + U(NPDE,NPTS), UE(NPDE,NPTS), + UOUT(NPDE,INTPTS,ITYPE), X(NPTS), XFIX(1), XI(1), + XOUT(INTPTS) INTEGER ISAVE(LISAVE) * .. External Functions .. DOUBLE PRECISION X01AAF EXTERNAL X01AAF * .. External Subroutines .. EXTERNAL BNDARY, D03PEK, D03PRF, D03PZF, EXACT, MONITF, + PDEDEF, UVINIT * .. Common blocks .. COMMON /PI/P * .. Executable Statements .. WRITE (NOUT,*) 'D03PRF Example Program Results' P = X01AAF(XX) ITRACE = 0 ITOL = 1 ATOL(1) = 0.5D-4 RTOL(1) = ATOL(1) * * Set remesh parameters * REMESH = .TRUE. NRMESH = 3 DXMESH = 0.0D0 CON = 5.0D0/(NPTS-1.0D0) XRATIO = 1.2D0 IPMINF = 0 * * Initialise mesh * DO 20 I = 1, NPTS X(I) = (I-1.0D0)/(NPTS-1.0D0) 20 CONTINUE * XOUT(1) = 0.0D0 XOUT(2) = 0.25D0 XOUT(3) = 0.5D0 XOUT(4) = 0.75D0 XOUT(5) = 1.0D0 * XI(1) = 0.0D0 NORM = 'A' LAOPT = 'F' IND = 0 ITASK = 1 * * Set THETA to .TRUE. if the Theta integrator is required * THETA = .FALSE. DO 40 I = 1, 30 ALGOPT(I) = 0.0D0 40 CONTINUE IF (THETA) THEN ALGOPT(1) = 2.0D0 ALGOPT(6) = 2.0D0 ALGOPT(7) = 1.0D0 END IF * * Loop over output value of t * TS = 0.0D0 TOUT = 0.0D0 * DO 60 IT = 1, 5 TOUT = 0.05D0*IT IFAIL = 1 * CALL D03PRF(NPDE,TS,TOUT,PDEDEF,BNDARY,UVINIT,U,NPTS,X,NLEFT, + NCODE,D03PEK,NXI,XI,NEQN,RTOL,ATOL,ITOL,NORM,LAOPT, + ALGOPT,REMESH,NXFIX,XFIX,NRMESH,DXMESH,TRMESH, + IPMINF,XRATIO,CON,MONITF,RSAVE,LRSAVE,ISAVE,LISAVE, + ITASK,ITRACE,IND,IFAIL) IF (IFAIL.NE.0) THEN WRITE (NOUT,99990) IFAIL GO TO 80 ELSE IF (IT.EQ.1) THEN WRITE (NOUT,99996) ATOL, NPTS WRITE (NOUT,99999) NRMESH WRITE (NOUT,99998) (XOUT(I),I=1,INTPTS) END IF * * Interpolate at output points .. CALL D03PZF(NPDE,0,U,NPTS,X,XOUT,INTPTS,ITYPE,UOUT,IFAIL) * Check against exact solution .. CALL EXACT(TS,NPDE,INTPTS,XOUT,UE) * WRITE (NOUT,99997) TS WRITE (NOUT,99994) (UOUT(1,I,1),I=1,INTPTS) WRITE (NOUT,99993) (UE(1,I),I=1,INTPTS) WRITE (NOUT,99992) (UOUT(2,I,1),I=1,INTPTS) WRITE (NOUT,99991) (UE(2,I),I=1,INTPTS) * 60 CONTINUE WRITE (NOUT,99995) ISAVE(1), ISAVE(2), ISAVE(3), ISAVE(5) 80 CONTINUE * 99999 FORMAT (' Remeshing every ',I3,' time steps',/) 99998 FORMAT (' X ',5F10.4,/) 99997 FORMAT (' T = ',F6.3) 99996 FORMAT (//' Accuracy requirement =',E10.3,' Number of points = ', + I3,/) 99995 FORMAT (' Number of integration steps in time = ',I6,/' Number o', + 'f function evaluations = ',I6,/' Number of Jacobian eval', + 'uations =',I6,/' Number of iterations = ',I6) 99994 FORMAT (' Approx U1',5F10.4) 99993 FORMAT (' Exact U1',5F10.4) 99992 FORMAT (' Approx U2',5F10.4) 99991 FORMAT (' Exact U2',5F10.4,/) 99990 FORMAT (1X,/1X,' ** D03PRF returned with IFAIL = ',I5) END * SUBROUTINE UVINIT(NPDE,NPTS,NXI,X,XI,U,NCODE,V) * .. Scalar Arguments .. INTEGER NCODE, NPDE, NPTS, NXI * .. Array Arguments .. DOUBLE PRECISION U(NPDE,NPTS), V(NCODE), X(NPTS), XI(NXI) * .. Scalars in Common .. DOUBLE PRECISION P * .. Local Scalars .. INTEGER I * .. Intrinsic Functions .. INTRINSIC EXP, SIN * .. Common blocks .. COMMON /PI/P * .. Executable Statements .. DO 20 I = 1, NPTS U(1,I) = EXP(X(I)) U(2,I) = X(I)**2 + SIN(2.0D0*P*X(I)**2) 20 CONTINUE RETURN END * SUBROUTINE PDEDEF(NPDE,T,X,U,UDOT,DUDX,NCODE,V,VDOT,RES,IRES) * .. Scalar Arguments .. DOUBLE PRECISION T, X INTEGER IRES, NCODE, NPDE * .. Array Arguments .. DOUBLE PRECISION DUDX(NPDE), RES(NPDE), U(NPDE), UDOT(NPDE), + V(NCODE), VDOT(NCODE) * .. Executable Statements .. IF (IRES.EQ.-1) THEN RES(1) = UDOT(1) RES(2) = UDOT(2) ELSE RES(1) = UDOT(1) + DUDX(1) + DUDX(2) RES(2) = UDOT(2) + 4.0D0*DUDX(1) + DUDX(2) END IF RETURN END * SUBROUTINE BNDARY(NPDE,T,IBND,NOBC,U,UDOT,NCODE,V,VDOT,RES,IRES) * .. Scalar Arguments .. DOUBLE PRECISION T INTEGER IBND, IRES, NCODE, NOBC, NPDE * .. Array Arguments .. DOUBLE PRECISION RES(NOBC), U(NPDE), UDOT(NPDE), V(NCODE), + VDOT(NCODE) * .. Scalars in Common .. DOUBLE PRECISION P * .. Local Scalars .. DOUBLE PRECISION PP * .. Intrinsic Functions .. INTRINSIC EXP, SIN * .. Common blocks .. COMMON /PI/P * .. Executable Statements .. PP = 2.0D0*P IF (IBND.EQ.0) THEN IF (IRES.EQ.-1) THEN RES(1) = 0.0D0 ELSE RES(1) = U(1) - 0.5D0*(EXP(T)+EXP(-3.0D0*T)) - + 0.25D0*(SIN(PP*9.0D0*T**2)-SIN(PP*T**2)) - + 2.0D0*T**2 END IF ELSE IF (IRES.EQ.-1) THEN RES(1) = 0.0D0 ELSE RES(1) = U(2) - (EXP(1.0D0-3.0D0*T)-EXP(1.0D0+T) + +0.5D0*(SIN(PP*(1.0D0-3.0D0*T)**2)+SIN(PP*(1.0D0+T) + **2))+1.0D0+5.0D0*T**2-2.0D0*T) END IF END IF RETURN END * SUBROUTINE EXACT(T,NPDE,NPTS,X,U) * Exact solution (for comparison purposes) * .. Scalar Arguments .. DOUBLE PRECISION T INTEGER NPDE, NPTS * .. Array Arguments .. DOUBLE PRECISION U(NPDE,NPTS), X(NPTS) * .. Scalars in Common .. DOUBLE PRECISION P * .. Local Scalars .. DOUBLE PRECISION PP INTEGER I * .. Intrinsic Functions .. INTRINSIC EXP, SIN * .. Common blocks .. COMMON /PI/P * .. Executable Statements .. PP = 2.0D0*P DO 20 I = 1, NPTS U(1,I) = 0.5D0*(EXP(X(I)+T)+EXP(X(I)-3.0D0*T)) + + 0.25D0*(SIN(PP*(X(I)-3.0D0*T)**2)-SIN(PP*(X(I)+T)**2)) + + 2.0D0*T**2 - 2.0D0*X(I)*T U(2,I) = EXP(X(I)-3.0D0*T) - EXP(X(I)+T) + 0.5D0*(SIN(PP*(X(I) + -3.0D0*T)**2)+SIN(PP*(X(I)+T)**2)) + X(I)**2 + + 5.0D0*T**2 - 2.0D0*X(I)*T 20 CONTINUE RETURN END * SUBROUTINE MONITF(T,NPTS,NPDE,X,U,FMON) * .. Scalar Arguments .. DOUBLE PRECISION T INTEGER NPDE, NPTS * .. Array Arguments .. DOUBLE PRECISION FMON(NPTS), U(NPDE,NPTS), X(NPTS) * .. Local Scalars .. DOUBLE PRECISION D2X1, D2X2, H1, H2, H3 INTEGER I * .. Intrinsic Functions .. INTRINSIC ABS, MAX * .. Executable Statements .. DO 20 I = 2, NPTS - 1 H1 = X(I) - X(I-1) H2 = X(I+1) - X(I) H3 = 0.5D0*(X(I+1)-X(I-1)) * Second derivatives .. D2X1 = ABS(((U(1,I+1)-U(1,I))/H2-(U(1,I)-U(1,I-1))/H1)/H3) D2X2 = ABS(((U(2,I+1)-U(2,I))/H2-(U(2,I)-U(2,I-1))/H1)/H3) FMON(I) = MAX(D2X1,D2X2) 20 CONTINUE FMON(1) = FMON(2) FMON(NPTS) = FMON(NPTS-1) RETURN END