! D02TXF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02txfe_mod ! D02TXF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: half = 0.5_nag_wp REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: three = 3.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: two = 2.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: xsplit = 30.0_nag_wp INTEGER, PARAMETER :: m1 = 3, m2 = 2, mmax = 3, & neq = 2, nin = 5, nlbc = 3, & nleft = 15, nout = 6, nrbc = 2, & nright = 10 ! .. Local Scalars .. REAL (KIND=nag_wp) :: el, en, s CONTAINS SUBROUTINE ffun(x,y,neq,m,f) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: f(neq) REAL (KIND=nag_wp), INTENT (IN) :: y(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Local Scalars .. REAL (KIND=nag_wp) :: t1, y11, y20 ! .. Executable Statements .. t1 = half*(three-en)*y(1,0) y11 = y(1,1) y20 = y(2,0) f(1) = (el**3)*(one-y20**2) + (el**2)*s*y11 - & el*(t1*y(1,2)+en*y11**2) f(2) = (el**2)*s*(y20-one) - el*(t1*y(2,1)+(en-one)*y11*y20) RETURN END SUBROUTINE ffun SUBROUTINE fjac(x,y,neq,m,dfdy) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: dfdy(neq,neq,0:*) REAL (KIND=nag_wp), INTENT (IN) :: y(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. dfdy(1,2,0) = -two*el**3*y(2,0) dfdy(1,1,0) = -el*half*(three-en)*y(1,2) dfdy(1,1,1) = el**2*s - el*two*en*y(1,1) dfdy(1,1,2) = -el*half*(three-en)*y(1,0) dfdy(2,2,0) = el**2*s - el*(en-one)*y(1,1) dfdy(2,2,1) = -el*half*(three-en)*y(1,0) dfdy(2,1,0) = -el*half*(three-en)*y(2,1) dfdy(2,1,1) = -el*(en-one)*y(2,0) RETURN END SUBROUTINE fjac SUBROUTINE gafun(ya,neq,m,nlbc,ga) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: neq, nlbc ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: ga(nlbc) REAL (KIND=nag_wp), INTENT (IN) :: ya(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. ga(1) = ya(1,0) ga(2) = ya(1,1) ga(3) = ya(2,0) RETURN END SUBROUTINE gafun SUBROUTINE gbfun(yb,neq,m,nrbc,gb) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: neq, nrbc ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: gb(nrbc) REAL (KIND=nag_wp), INTENT (IN) :: yb(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. gb(1) = yb(1,1) gb(2) = yb(2,0) - one RETURN END SUBROUTINE gbfun SUBROUTINE gajac(ya,neq,m,nlbc,dgady) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: neq, nlbc ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: dgady(nlbc,neq,0:*) REAL (KIND=nag_wp), INTENT (IN) :: ya(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. dgady(1,1,0) = one dgady(2,1,1) = one dgady(3,2,0) = one RETURN END SUBROUTINE gajac SUBROUTINE gbjac(yb,neq,m,nrbc,dgbdy) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: neq, nrbc ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: dgbdy(nrbc,neq,0:*) REAL (KIND=nag_wp), INTENT (IN) :: yb(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Executable Statements .. dgbdy(1,1,1) = one dgbdy(2,2,0) = one RETURN END SUBROUTINE gbjac SUBROUTINE guess(x,neq,m,y,dym) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: dym(neq) REAL (KIND=nag_wp), INTENT (INOUT) :: y(neq,0:*) INTEGER, INTENT (IN) :: m(neq) ! .. Local Scalars .. REAL (KIND=nag_wp) :: ex, expmx ! .. Intrinsic Functions .. INTRINSIC exp ! .. Executable Statements .. ex = x*el expmx = exp(-ex) y(1,0) = -ex**2*expmx y(1,1) = (-two*ex+ex**2)*expmx y(1,2) = (-two+4.0E0_nag_wp*ex-ex**2)*expmx y(2,0) = one - expmx y(2,1) = expmx dym(1) = (6.0E0_nag_wp-6.0E0_nag_wp*ex+ex**2)*expmx dym(2) = -expmx RETURN END SUBROUTINE guess END MODULE d02txfe_mod PROGRAM d02txfe ! D02TXF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d02tkf, d02tvf, d02txf, d02tyf, d02tzf USE d02txfe_mod, ONLY : el, en, ffun, fjac, gafun, gajac, gbfun, gbjac, & guess, m1, m2, mmax, nag_wp, neq, nin, nlbc, & nleft, nout, nrbc, nright, one, s, two, xsplit ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: dx, el_init, ermx, s_init, xx INTEGER :: i, iermx, ifail, ijermx, j, & liwork, lrwork, mxmesh, ncol, & ncont, nmesh ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: mesh(:), rwork(:) REAL (KIND=nag_wp) :: tol(neq), y(neq,0:mmax-1) INTEGER, ALLOCATABLE :: ipmesh(:), iwork(:) INTEGER :: m(neq) ! .. Intrinsic Functions .. INTRINSIC min, real ! .. Executable Statements .. WRITE (nout,*) 'D02TXF Example Program Results' WRITE (nout,*) ! Skip heading in data file READ (nin,*) ! Read method parameters READ (nin,*) ncol, nmesh, mxmesh READ (nin,*) tol(1:neq) liwork = mxmesh*(11*neq+6) lrwork = mxmesh*(109*neq**2+78*neq+7) ALLOCATE (mesh(mxmesh),rwork(lrwork),ipmesh(mxmesh),iwork(liwork)) ! Read problem (initial) parameters READ (nin,*) en, el_init, s_init ! Initialise data el = el_init s = s_init m(1) = m1 m(2) = m2 dx = one/real(nmesh-1,kind=nag_wp) mesh(1) = 0.0_nag_wp DO i = 2, nmesh - 1 mesh(i) = mesh(i-1) + dx END DO mesh(nmesh) = 1.0_nag_wp ipmesh(1) = 1 ipmesh(2:nmesh-1) = 2 ipmesh(nmesh) = 1 ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 CALL d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,rwork, & lrwork,iwork,liwork,ifail) ! Initialize number of continuation steps in el and s READ (nin,*) ncont CONT: DO j = 1, ncont WRITE (nout,99997) tol(1), el, s ! Solve ifail = -1 CALL d02tkf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rwork,iwork, & ifail) IF (ifail/=0) EXIT CONT ! Extract mesh ifail = 0 CALL d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rwork,iwork, & ifail) WRITE (nout,99996) nmesh, ermx, iermx, ijermx ! Print solution components on mesh WRITE (nout,99999) ! Left side domain [0,xsplit], evaluate at nleft+1 uniform grid points. dx = xsplit/real(nleft,kind=nag_wp)/el xx = 0.0_nag_wp DO i = 0, nleft ifail = 0 CALL d02tyf(xx,y,neq,mmax,rwork,iwork,ifail) WRITE (nout,99998) xx*el, y(1,0), y(2,0) xx = xx + dx END DO ! Right side domain (xsplit,L], evaluate at nright uniform grid points. dx = (el-xsplit)/real(nright,kind=nag_wp)/el xx = xsplit/el DO i = 1, nright xx = min(1.0_nag_wp,xx+dx) ifail = 0 CALL d02tyf(xx,y,neq,mmax,rwork,iwork,ifail) WRITE (nout,99998) xx*el, y(1,0), y(2,0) END DO ! Select mesh for continuation and update continuation parameters. IF (j