! D02TVF Example Program Text
! Mark 26.2 Release. NAG Copyright 2017.
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
! .. Accessibility Statements ..
Private
Public :: ffun, fjac, gafun, gajac, gbfun, &
gbjac, guess
! .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: 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, Public :: mmax = 1, neq = 6, nin = 5, &
nlbc = 3, nout = 6, nrbc = 3
! .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: beta0, eta, lambda, mu
! .. Local Arrays ..
Integer, Public, Save :: m(neq) = (/1,1,1,1,1,1/)
Contains
Subroutine ffun(x,y,neq,m,f,iuser,ruser)
! .. Use Statements ..
Use nag_library, Only: x01aaf
! .. 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 (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: y(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
Integer, Intent (In) :: m(neq)
! .. Local Scalars ..
Real (Kind=nag_wp) :: beta
! .. Intrinsic Procedures ..
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,iuser,ruser)
! .. Use Statements ..
Use nag_library, Only: x01aaf
! .. 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:*), ruser(*)
Real (Kind=nag_wp), Intent (In) :: y(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
Integer, Intent (In) :: m(neq)
! .. Local Scalars ..
Real (Kind=nag_wp) :: beta
! .. Intrinsic Procedures ..
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,iuser,ruser)
! .. Scalar Arguments ..
Integer, Intent (In) :: neq, nlbc
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: ga(nlbc)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: ya(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
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,iuser,ruser)
! .. Scalar Arguments ..
Integer, Intent (In) :: neq, nrbc
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: gb(nrbc)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: yb(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
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,iuser,ruser)
! .. Scalar Arguments ..
Integer, Intent (In) :: neq, nlbc
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: dgady(nlbc,neq,0:*), ruser(*)
Real (Kind=nag_wp), Intent (In) :: ya(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
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,iuser,ruser)
! .. Scalar Arguments ..
Integer, Intent (In) :: neq, nrbc
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: dgbdy(nrbc,neq,0:*), ruser(*)
Real (Kind=nag_wp), Intent (In) :: yb(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
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,iuser,ruser)
! .. 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) :: ruser(*), y(neq,0:*)
Integer, Intent (Inout) :: iuser(*)
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 d02tvfe_mod, Only: beta0, eta, ffun, fjac, gafun, gajac, gbfun, &
gbjac, guess, lambda, m, mmax, mu, neq, nin, &
nlbc, nout, nrbc, one
Use nag_library, Only: d02tlf, d02tvf, d02tyf, d02tzf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: dx, ermx
Integer :: i, iermx, ifail, ijermx, licomm, &
lrcomm, mxmesh, ncol, nmesh
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: mesh(:), rcomm(:), tols(:), y(:,:)
Real (Kind=nag_wp) :: ruser(1)
Integer, Allocatable :: icomm(:), ipmesh(:)
Integer :: iuser(2)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*) 'D02TVF Example Program Results'
Write (nout,*)
! Skip heading in data file
Read (nin,*)
Read (nin,*) ncol, nmesh, mxmesh
Allocate (mesh(mxmesh),tols(neq),y(neq,0:mmax-1),ipmesh(mxmesh))
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
! Workspace query to get size of rcomm and icomm
ifail = 0
Call d02tvf(neq,m,nlbc,nrbc,ncol,tols,mxmesh,nmesh,mesh,ipmesh,ruser,0, &
iuser,2,ifail)
lrcomm = iuser(1)
licomm = iuser(2)
Allocate (rcomm(lrcomm),icomm(licomm))
! Initialize
ifail = 0
Call d02tvf(neq,m,nlbc,nrbc,ncol,tols,mxmesh,nmesh,mesh,ipmesh,rcomm, &
lrcomm,icomm,licomm,ifail)
! Solve
ifail = -1
Call d02tlf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rcomm,icomm,iuser, &
ruser,ifail)
! Extract mesh.
ifail = -1
Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rcomm,icomm, &
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,rcomm,icomm,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