! D02TKF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02tkfe_mod ! D02TKF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: mmax = 3, neq = 3, nin = 5, & nlbc = 3, nout = 6, nrbc = 3 ! .. Local Scalars .. REAL (KIND=nag_wp) :: omega REAL (KIND=nag_wp) :: one = 1.0_nag_wp REAL (KIND=nag_wp) :: sqrofr ! .. Local Arrays .. INTEGER :: m(neq) = (/ 1, 3, 2/) 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) ! .. Executable Statements .. f(1) = y(2,0) f(2) = -(y(1,0)*y(2,2)+y(3,0)*y(3,1))*sqrofr f(3) = (y(2,0)*y(3,0)-y(1,0)*y(3,1))*sqrofr 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) = one dfdy(2,1,0) = -y(2,2)*sqrofr dfdy(2,2,2) = -y(1,0)*sqrofr dfdy(2,3,0) = -y(3,1)*sqrofr dfdy(2,3,1) = -y(3,0)*sqrofr dfdy(3,1,0) = -y(3,1)*sqrofr dfdy(3,2,0) = y(3,0)*sqrofr dfdy(3,3,0) = y(2,0)*sqrofr dfdy(3,3,1) = -y(1,0)*sqrofr 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(2,0) ga(3) = ya(3,0) - omega 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) gb(2) = yb(2,0) gb(3) = yb(3,0) + omega 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,2,0) = one dgady(3,3,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(2,2,0) = one dgbdy(3,3,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) :: five = 5.0_nag_wp REAL (KIND=nag_wp) :: half = 0.5_nag_wp REAL (KIND=nag_wp) :: ten = 10.0_nag_wp REAL (KIND=nag_wp) :: two = 2.0_nag_wp ! .. Executable Statements .. y(1,0) = -(x-half)*(x*(x-one))**2 y(2,0) = -x*(x-one)*(five*x*(x-one)+one) y(2,1) = -(two*x-one)*(ten*x*(x-one)+one) y(2,2) = -12.0_nag_wp*(five*x*(x-one)+x) y(3,0) = -8.0_nag_wp*omega*(x-half)**3 y(3,1) = -24.0_nag_wp*omega*(x-half)**2 dym(1) = y(2,0) dym(2) = -120.0_nag_wp*(x-half) dym(3) = -56.0_nag_wp*omega*(x-half) RETURN END SUBROUTINE guess END MODULE d02tkfe_mod PROGRAM d02tkfe ! D02TKF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d02tkf, d02tvf, d02txf, d02tyf, d02tzf USE d02tkfe_mod, ONLY : ffun, fjac, gafun, gajac, gbfun, gbjac, guess, & m, mmax, nag_wp, neq, nin, nlbc, nout, nrbc, & omega, one, sqrofr ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: dx, ermx, r INTEGER :: i, iermx, ifail, ijermx, j, & liwork, lrwork, mxmesh, ncol, & ncont, nmesh ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: mesh(:), tol(:), work(:), y(:,:) INTEGER, ALLOCATABLE :: ipmesh(:), iwork(:) ! .. Intrinsic Functions .. INTRINSIC real, sqrt ! .. Executable Statements .. WRITE (nout,*) 'D02TKF 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),tol(neq),work(lrwork),y(neq,0:mmax-1), & ipmesh(mxmesh),iwork(liwork)) READ (nin,*) omega READ (nin,*) tol(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 ! Initial integrator for given problem. ifail = 0 CALL d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,work, & lrwork,iwork,liwork,ifail) ! Number of continuation steps (last r=100**ncont, sqrofr=10**ncont) READ (nin,*) ncont ! Initialize problem continuation parameter. READ (nin,*) r sqrofr = sqrt(r) CONTN: DO j = 1, ncont WRITE (nout,99999) tol(1), r ! Solve problem. ifail = -1 CALL d02tkf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,work,iwork, & ifail) ! Extract mesh ifail = -1 CALL d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,work,iwork, & ifail) IF (ifail==1) THEN EXIT CONTN END IF ! Print mesh and error statistics. WRITE (nout,99998) nmesh, ermx, iermx, ijermx WRITE (nout,99997) (i,ipmesh(i),mesh(i),i=1,nmesh) ! Print solution components on mesh. WRITE (nout,99996) DO i = 1, nmesh ifail = 0 CALL d02tyf(mesh(i),y,neq,mmax,work,iwork,ifail) WRITE (nout,99995) mesh(i), y(1:neq,0) END DO IF (j==ncont) EXIT CONTN ! Modify continuation parameter. r = 100.0_nag_wp*r sqrofr = sqrt(r) ! Select mesh for continuation and call continuation primer routine. ipmesh(2:nmesh-1) = 2 ifail = 0 CALL d02txf(mxmesh,nmesh,mesh,ipmesh,work,iwork,ifail) END DO CONTN 99999 FORMAT (/' Tolerance = ',1P,E8.1,' R = ',E10.3) 99998 FORMAT (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, & ' in interval ',I4,' for component ',I4/) 99997 FORMAT (/' Mesh points:'/4(I4,'(',I1,')',E11.4)) 99996 FORMAT (/' x f f'' g') 99995 FORMAT (' ',F8.3,1X,3F9.4) END PROGRAM d02tkfe