! D02TKF Example Program Text
! Mark 25 Release. NAG Copyright 2014.
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
! .. Accessibility Statements ..
Private
Public :: ffun, fjac, gafun, gajac, gbfun, &
gbjac, guess
! .. Parameters ..
Integer, Parameter, Public :: mmax = 3, neq = 3, nin = 5, &
nlbc = 3, nout = 6, nrbc = 3
! .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: omega
Real (Kind=nag_wp), Public, Save :: one = 1.0_nag_wp
Real (Kind=nag_wp), Public, Save :: sqrofr
! .. Local Arrays ..
Integer, Public, Save :: m(neq) = (/1,3,2/)
Contains
Subroutine ffun(x,y,neq,m,f)
! .. 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)
! .. 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)
! .. 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)
! .. 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)
! .. 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)
! .. 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)
! .. 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,0) = -(x-0.5_nag_wp)*(x*(x-one))**2
y(2,0) = -x*(x-one)*(5._nag_wp*x*(x-one)+one)
y(2,1) = -(2._nag_wp*x-one)*(10._nag_wp*x*(x-one)+one)
y(2,2) = -12.0_nag_wp*(5._nag_wp*x*(x-one)+x)
y(3,0) = -8.0_nag_wp*omega*(x-0.5_nag_wp)**3
y(3,1) = -24.0_nag_wp*omega*(x-0.5_nag_wp)**2
dym(1) = y(2,0)
dym(2) = -120.0_nag_wp*(x-0.5_nag_wp)
dym(3) = -56.0_nag_wp*omega*(x-0.5_nag_wp)
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, nag_wp
Use d02tkfe_mod, Only: ffun, fjac, gafun, gajac, gbfun, gbjac, guess, m, &
mmax, 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, &
licomm, lrcomm, mxmesh, ncol, &
ncont, nmesh
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: mesh(:), tol(:), work(:), y(:,:)
Integer, Allocatable :: icomm(:), ipmesh(:)
! .. Intrinsic Procedures ..
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
licomm = 23 + neq + mxmesh
lrcomm = mxmesh*(109*neq**2+78*neq+7)
Allocate (mesh(mxmesh),tol(neq),work(lrcomm),y(neq,0:mmax-1), &
ipmesh(mxmesh),icomm(licomm))
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, &
lrcomm,icomm,licomm,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,icomm,ifail)
! Extract mesh
ifail = -1
Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,work,icomm, &
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,icomm,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,icomm,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