* D03PJA Example Program Text * Mark 20 Release. NAG Copyright 2001. * .. Parameters .. INTEGER NOUT PARAMETER (NOUT=6) INTEGER NBKPTS, NEL, NPDE, NPOLY, NPTS, NCODE, M, NXI, + NEQN, NIW, NPL1, NWKRES, LENODE, NW PARAMETER (NBKPTS=11,NEL=NBKPTS-1,NPDE=1,NPOLY=2, + NPTS=NEL*NPOLY+1,NCODE=1,M=0,NXI=1, + NEQN=NPDE*NPTS+NCODE,NIW=24,NPL1=NPOLY+1, + NWKRES=3*NPL1*NPL1+NPL1* + (NPDE*NPDE+6*NPDE+NBKPTS+1)+8*NPDE+NXI*(5*NPDE+1) + +NCODE+3,LENODE=11*NEQN+50, + NW=NEQN*NEQN+NEQN+NWKRES+LENODE) * .. 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(NBKPTS), RTOL(1), + RUSER(1), RWSAV(1100), U(NEQN), W(NW), X(NPTS), + XBKPTS(NBKPTS), XI(1) INTEGER IUSER(1), IW(NIW), IWSAV(505) LOGICAL LWSAV(100) CHARACTER*80 CWSAV(10) * .. External Subroutines .. EXTERNAL BNDARY, D03PJA, EXACT, ODEDEF, PDEDEF, UVINIT * .. Executable Statements .. WRITE (NOUT,*) 'D03PJA Example Program Results' ITRACE = 0 ITOL = 1 ATOL(1) = 1.0D-4 RTOL(1) = ATOL(1) * * Set break-points * DO 20 I = 1, NBKPTS XBKPTS(I) = (I-1.0D0)/(NBKPTS-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 * * Loop over output value of t * RUSER(1) = 1.0D-4 TOUT = 0.0D0 DO 60 IT = 1, 5 TOUT = 0.1D0*(2**IT) IFAIL = 1 * CALL D03PJA(NPDE,M,RUSER(1),TOUT,PDEDEF,BNDARY,U,NBKPTS,XBKPTS, + NPOLY,NPTS,X,NCODE,ODEDEF,NXI,XI,NEQN,UVINIT,RTOL, + ATOL,ITOL,NORM,LAOPT,ALGOPT,W,NW,IW,NIW,ITASK, + ITRACE,IND,IUSER,RUSER,CWSAV,LWSAV,IWSAV,RWSAV, + IFAIL) * IF (IFAIL.NE.0) THEN WRITE (NOUT,99992) IFAIL GO TO 80 ELSE IF (IT.EQ.1) THEN WRITE (NOUT,99999) NPOLY, NEL WRITE (NOUT,99996) ATOL, NPTS WRITE (NOUT,99998) (XBKPTS(I),I=1,7,2), XBKPTS(11) END IF * * Check against the exact solution * CALL EXACT(TOUT,NBKPTS,XBKPTS,EXY) WRITE (NOUT,99997) RUSER(1) WRITE (NOUT,99994) (U(I),I=1,13,4), U(21), U(22) WRITE (NOUT,99993) (EXY(I),I=1,7,2), EXY(11), RUSER(1) 60 CONTINUE WRITE (NOUT,99995) IW(1), IW(2), IW(3), IW(5) * 80 CONTINUE * 99999 FORMAT (' Degree of Polynomial =',I4,' No. of elements =',I4,/) 99998 FORMAT (' X ',5F9.3,/) 99997 FORMAT (' T = ',F6.3) 99996 FORMAT (//' Simple coupled PDE using BDF ',/' Accuracy require', + 'ment =',E10.3,' Number of points = ',I4,/) 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 (1X,'App. sol. ',F7.3,4F9.3,' ODE sol. =',F8.3) 99993 FORMAT (1X,'Exact sol. ',F7.3,4F9.3,' ODE sol. =',F8.3,/) 99992 FORMAT (1X,/1X,' ** D03PJA returned with IFAIL = ',I5) END * SUBROUTINE UVINIT(NPDE,NPTS,X,U,NCODE,V,IUSER,RUSER) * Routine for PDE initial values (start time is 0.1D-6) * .. Scalar Arguments .. INTEGER NCODE, NPDE, NPTS * .. Array Arguments .. DOUBLE PRECISION RUSER(1), U(NPDE,NPTS), V(*), X(NPTS) INTEGER IUSER(1) * .. Local Scalars .. DOUBLE PRECISION TS INTEGER I * .. Intrinsic Functions .. INTRINSIC EXP * .. Executable Statements .. TS = RUSER(1) V(1) = TS DO 20 I = 1, NPTS U(1,I) = EXP(TS*(1.0D0-X(I))) - 1.0D0 20 CONTINUE RETURN END * SUBROUTINE ODEDEF(NPDE,T,NCODE,V,VDOT,NXI,XI,UCP,UCPX,RCP,UCPT, + UCPTX,F,IRES,IUSER,RUSER) * .. Scalar Arguments .. DOUBLE PRECISION T INTEGER IRES, NCODE, NPDE, NXI * .. Array Arguments .. DOUBLE PRECISION F(*), RCP(NPDE,*), RUSER(*), UCP(NPDE,*), + UCPT(NPDE,*), UCPTX(NPDE,*), UCPX(NPDE,*), V(*), + VDOT(*), XI(*) INTEGER IUSER(*) * .. Executable Statements .. IF (IRES.EQ.1) THEN F(1) = VDOT(1) - V(1)*UCP(1,1) - UCPX(1,1) - 1.0D0 - T ELSE IF (IRES.EQ.-1) THEN F(1) = VDOT(1) END IF RETURN END * SUBROUTINE PDEDEF(NPDE,T,X,NPTL,U,DUDX,NCODE,V,VDOT,P,Q,R,IRES, + IUSER,RUSER) * .. Scalar Arguments .. DOUBLE PRECISION T INTEGER IRES, NCODE, NPDE, NPTL * .. Array Arguments .. DOUBLE PRECISION DUDX(NPDE,NPTL), P(NPDE,NPDE,NPTL), Q(NPDE,NPTL), + R(NPDE,NPTL), RUSER(*), U(NPDE,NPTL), V(*), + VDOT(*), X(NPTL) INTEGER IUSER(*) * .. Local Scalars .. INTEGER I * .. Executable Statements .. DO 20 I = 1, NPTL P(1,1,I) = V(1)*V(1) R(1,I) = DUDX(1,I) Q(1,I) = -X(I)*DUDX(1,I)*V(1)*VDOT(1) 20 CONTINUE RETURN END * SUBROUTINE BNDARY(NPDE,T,U,UX,NCODE,V,VDOT,IBND,BETA,GAMMA,IRES, + IUSER,RUSER) * .. Scalar Arguments .. DOUBLE PRECISION T INTEGER IBND, IRES, NCODE, NPDE * .. Array Arguments .. DOUBLE PRECISION BETA(NPDE), GAMMA(NPDE), RUSER(*), U(NPDE), + UX(NPDE), V(*), VDOT(*) INTEGER IUSER(*) * .. Intrinsic Functions .. INTRINSIC EXP * .. Executable Statements .. BETA(1) = 1.0D0 IF (IBND.EQ.0) THEN GAMMA(1) = -V(1)*EXP(T) ELSE GAMMA(1) = -V(1)*VDOT(1) END IF RETURN END * SUBROUTINE EXACT(TIME,NPTS,X,U) * Exact solution (for comparison purposes) * .. Scalar Arguments .. DOUBLE PRECISION TIME INTEGER NPTS * .. Array Arguments .. DOUBLE PRECISION U(NPTS), X(NPTS) * .. Local Scalars .. INTEGER I * .. Intrinsic Functions .. INTRINSIC EXP * .. Executable Statements .. DO 20 I = 1, NPTS U(I) = EXP(TIME*(1.0D0-X(I))) - 1.0D0 20 CONTINUE RETURN END