! D02TVF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02tvfe_mod ! D02TVF 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 :: one = 1.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: two = 2.0_nag_wp REAL (KIND=nag_wp), PARAMETER :: zero = 0.0_nag_wp INTEGER, PARAMETER :: mmax = 1, neq = 6, nin = 5, & nlbc = 3, nout = 6, nrbc = 3 ! .. Local Scalars .. REAL (KIND=nag_wp) :: beta0, eta, lambda, mu ! .. Local Arrays .. INTEGER :: m(neq) = (/ 1, 1, 1, 1, 1, 1/) CONTAINS SUBROUTINE ffun(x,y,neq,m,f) ! .. Use Statements .. USE nag_library, ONLY : x01aaf ! .. 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) :: beta ! .. Intrinsic Functions .. INTRINSIC cos ! .. Executable Statements .. beta = beta0*(one+cos(two*x01aaf(beta)*x)) f(1) = mu - beta*y(1,0)*y(3,0) f(2) = beta*y(1,0)*y(3,0) - y(2,0)/lambda f(3) = y(2,0)/lambda - y(3,0)/eta f(4:6) = zero RETURN END SUBROUTINE ffun SUBROUTINE fjac(x,y,neq,m,dfdy) ! .. Use Statements .. USE nag_library, ONLY : x01aaf ! .. 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) ! .. Local Scalars .. REAL (KIND=nag_wp) :: beta ! .. Intrinsic Functions .. INTRINSIC cos ! .. Executable Statements .. beta = beta0*(one+cos(two*x01aaf(beta)*x)) dfdy(1,1,0) = -beta*y(3,0) dfdy(1,3,0) = -beta*y(1,0) dfdy(2,1,0) = beta*y(3,0) dfdy(2,2,0) = -one/lambda dfdy(2,3,0) = beta*y(1,0) dfdy(3,2,0) = one/lambda dfdy(3,3,0) = -one/eta 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) - ya(4,0) ga(2) = ya(2,0) - ya(5,0) ga(3) = ya(3,0) - ya(6,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,0) - yb(4,0) gb(2) = yb(2,0) - yb(5,0) gb(3) = yb(3,0) - yb(6,0) 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(1,4,0) = -one dgady(2,2,0) = one dgady(2,5,0) = -one dgady(3,3,0) = one dgady(3,6,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,0) = one dgbdy(1,4,0) = -one dgbdy(2,2,0) = one dgbdy(2,5,0) = -one dgbdy(3,3,0) = one dgbdy(3,6,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) ! .. Executable Statements .. y(1:3,0) = one y(4,0) = y(1,0) y(5,0) = y(2,0) y(6,0) = y(3,0) dym(1:neq) = zero RETURN END SUBROUTINE guess END MODULE d02tvfe_mod PROGRAM d02tvfe ! D02TVF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d02tkf, d02tvf, d02tyf, d02tzf USE d02tvfe_mod, ONLY : beta0, eta, ffun, fjac, gafun, gajac, gbfun, & gbjac, guess, lambda, m, mmax, mu, nag_wp, neq, & nin, nlbc, nout, nrbc, one ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: dx, ermx INTEGER :: i, iermx, ifail, ijermx, liwork, & lrwork, mxmesh, ncol, nmesh ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: mesh(:), rwork(:), tols(:), y(:,:) INTEGER, ALLOCATABLE :: ipmesh(:), iwork(:) ! .. Intrinsic Functions .. INTRINSIC real ! .. Executable Statements .. WRITE (nout,*) 'D02TVF Example Program Results' WRITE (nout,*) ! Skip heading in data file READ (nin,*) READ (nin,*) ncol, nmesh, mxmesh liwork = mxmesh*(11*neq+6) lrwork = mxmesh*(109*neq**2+78*neq+7) ALLOCATE (mesh(mxmesh),tols(neq),rwork(lrwork),y(neq,0:mmax-1), & ipmesh(mxmesh),iwork(liwork)) READ (nin,*) beta0, eta, lambda, mu READ (nin,*) tols(1:neq) 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) = one ipmesh(1) = 1 ipmesh(2:nmesh-1) = 2 ipmesh(nmesh) = 1 ! Initialize ifail = 0 CALL d02tvf(neq,m,nlbc,nrbc,ncol,tols,mxmesh,nmesh,mesh,ipmesh,rwork, & lrwork,iwork,liwork,ifail) ! Solve ifail = -1 CALL d02tkf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rwork,iwork,ifail) ! Extract mesh. ifail = -1 CALL d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rwork,iwork, & ifail) IF (ifail/=1) THEN ! Print mesh statistics WRITE (nout,99999) nmesh, ermx, iermx, ijermx WRITE (nout,99998) (i,ipmesh(i),mesh(i),i=1,nmesh) ! Print solution on mesh. WRITE (nout,99997) DO i = 1, nmesh ifail = 0 CALL d02tyf(mesh(i),y,neq,mmax,rwork,iwork,ifail) WRITE (nout,99996) mesh(i), y(1:3,0) END DO END IF 99999 FORMAT (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, & ' in interval ',I4,' for component ',I4/) 99998 FORMAT (/' Mesh points:'/4(I4,'(',I1,')',F7.4)) 99997 FORMAT (/' Computed solution at mesh points'/' x y1 ', & ' y2 y3') 99996 FORMAT (1X,F6.3,1X,3E11.3) END PROGRAM d02tvfe