* D03PKF Example Program Text * Mark 16 Release. NAG Copyright 1993. * .. Parameters .. INTEGER NOUT PARAMETER (NOUT=6) INTEGER NPDE, NPTS, NCODE, NXI, NLEFT, NEQN, LISAVE, + NWKRES, LENODE, LRSAVE PARAMETER (NPDE=2,NPTS=21,NCODE=1,NXI=1,NLEFT=1, + NEQN=NPDE*NPTS+NCODE,LISAVE=24, + NWKRES=NPDE*(NPTS+6*NXI+3*NPDE+15) + +NCODE+NXI+7*NPTS+2,LENODE=11*NEQN+50, + LRSAVE=NEQN*NEQN+NEQN+NWKRES+LENODE) * .. Scalars in Common .. DOUBLE PRECISION TS * .. Local Scalars .. DOUBLE PRECISION TOUT INTEGER I, IFAIL, IND, IT, ITASK, ITOL, ITRACE LOGICAL THETA CHARACTER LAOPT, NORM * .. Local Arrays .. DOUBLE PRECISION ALGOPT(30), ATOL(1), EXY(NEQN), RSAVE(LRSAVE), + RTOL(1), U(NEQN), X(NPTS), XI(1) INTEGER ISAVE(LISAVE) * .. External Subroutines .. EXTERNAL BNDARY, D03PKF, EXACT, ODEDEF, PDEDEF, UVINIT * .. Common blocks .. COMMON /TAXIS/TS * .. Executable Statements .. WRITE (NOUT,*) 'D03PKF Example Program Results' ITRACE = 0 ITOL = 1 ATOL(1) = 0.1D-3 RTOL(1) = ATOL(1) * * Set spatial-mesh points * DO 20 I = 1, NPTS X(I) = (I-1.0D0)/(NPTS-1.0D0) 20 CONTINUE * XI(1) = 1.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 ELSE ALGOPT(1) = 0.0D0 END IF ALGOPT(1) = 1.0D0 ALGOPT(13) = 0.5D-2 * * Loop over output value of t * TS = 1.0D-4 TOUT = 0.0D0 * CALL UVINIT(NPDE,NPTS,X,U,NCODE,NEQN) * DO 60 IT = 1, 5 TOUT = 0.1D0*(2**IT) IFAIL = 1 * CALL D03PKF(NPDE,TS,TOUT,PDEDEF,BNDARY,U,NPTS,X,NLEFT,NCODE, + ODEDEF,NXI,XI,NEQN,RTOL,ATOL,ITOL,NORM,LAOPT, + ALGOPT,RSAVE,LRSAVE,ISAVE,LISAVE,ITASK,ITRACE,IND, + IFAIL) * IF (IFAIL.NE.0) THEN WRITE (NOUT,99993) IFAIL GO TO 80 ELSE IF (IT.EQ.1) THEN WRITE (NOUT,99997) ATOL, NPTS WRITE (NOUT,99999) X(1), X(5), X(9), X(13), X(21) END IF * * Check against the exact solution * CALL EXACT(TOUT,NEQN,NPTS,X,EXY) * WRITE (NOUT,99998) TS WRITE (NOUT,99995) (U(I),I=1,25,8), U(41), U(43) WRITE (NOUT,99994) (EXY(I),I=1,25,8), EXY(41), TS 60 CONTINUE WRITE (NOUT,99996) ISAVE(1), ISAVE(2), ISAVE(3), ISAVE(5) * 80 CONTINUE * 99999 FORMAT (' X ',5F9.3,/) 99998 FORMAT (' T = ',F6.3) 99997 FORMAT (//' Accuracy requirement =',E10.3,' Number of points = ', + I3,/) 99996 FORMAT (' Number of integration steps in time = ',I6,/' Number o', + 'f function evaluations = ',I6,/' Number of Jacobian eval', + 'uations =',I6,/' Number of iterations = ',I6) 99995 FORMAT (1X,'App. sol. ',F7.3,4F9.3,' ODE sol. =',F8.3) 99994 FORMAT (1X,'Exact sol. ',F7.3,4F9.3,' ODE sol. =',F8.3,/) 99993 FORMAT (1X,/1X,' ** D03PKF returned with IFAIL = ',I5) END * SUBROUTINE UVINIT(NPDE,NPTS,X,U,NCODE,NEQN) * Routine for PDE initial values * .. Scalar Arguments .. INTEGER NCODE, NEQN, NPDE, NPTS * .. Array Arguments .. DOUBLE PRECISION U(NEQN), X(NPTS) * .. Scalars in Common .. DOUBLE PRECISION TS * .. Local Scalars .. INTEGER I, K * .. Intrinsic Functions .. INTRINSIC EXP * .. Common blocks .. COMMON /TAXIS/TS * .. Executable Statements .. K = 1 DO 20 I = 1, NPTS U(K) = EXP(TS*(1.0D0-X(I))) - 1.0D0 U(K+1) = -TS*EXP(TS*(1.0D0-X(I))) K = K + 2 20 CONTINUE U(NEQN) = TS RETURN END * SUBROUTINE ODEDEF(NPDE,T,NCODE,V,VDOT,NXI,XI,UCP,UCPX,UCPT,F,IRES) * .. Scalar Arguments .. DOUBLE PRECISION T INTEGER IRES, NCODE, NPDE, NXI * .. Array Arguments .. DOUBLE PRECISION F(*), UCP(NPDE,*), UCPT(NPDE,*), UCPX(NPDE,*), + V(*), VDOT(*), XI(*) * .. Executable Statements .. IF (IRES.EQ.-1) THEN F(1) = VDOT(1) ELSE F(1) = VDOT(1) - V(1)*UCP(1,1) - UCP(2,1) - 1.0D0 - T END IF RETURN END * SUBROUTINE PDEDEF(NPDE,T,X,U,UDOT,UX,NCODE,V,VDOT,RES,IRES) * .. Scalar Arguments .. DOUBLE PRECISION T, X INTEGER IRES, NCODE, NPDE * .. Array Arguments .. DOUBLE PRECISION RES(NPDE), U(NPDE), UDOT(NPDE), UX(NPDE), V(*), + VDOT(*) * .. Executable Statements .. IF (IRES.EQ.-1) THEN RES(1) = V(1)*V(1)*UDOT(1) - X*U(2)*V(1)*VDOT(1) RES(2) = 0.0D0 ELSE RES(1) = V(1)*V(1)*UDOT(1) - X*U(2)*V(1)*VDOT(1) - UX(2) RES(2) = U(2) - UX(1) 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(*), VDOT(*) * .. Intrinsic Functions .. INTRINSIC EXP * .. Executable Statements .. IF (IBND.EQ.0) THEN IF (IRES.EQ.-1) THEN RES(1) = 0.0D0 ELSE RES(1) = U(2) + V(1)*EXP(T) END IF ELSE IF (IRES.EQ.-1) THEN RES(1) = V(1)*VDOT(1) ELSE RES(1) = U(2) + V(1)*VDOT(1) END IF END IF RETURN END * SUBROUTINE EXACT(TIME,NEQN,NPTS,X,U) * Exact solution (for comparison purposes) * .. Scalar Arguments .. DOUBLE PRECISION TIME INTEGER NEQN, NPTS * .. Array Arguments .. DOUBLE PRECISION U(NEQN), X(NPTS) * .. Local Scalars .. INTEGER I, K * .. Intrinsic Functions .. INTRINSIC EXP * .. Executable Statements .. K = 1 DO 20 I = 1, NPTS U(K) = EXP(TIME*(1.0D0-X(I))) - 1.0D0 K = K + 2 20 CONTINUE RETURN END