! D02TXF Example Program Text
! Mark 27.3 Release. NAG Copyright 2021.
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
! .. Accessibility Statements ..
Private
Public :: ffun, fjac, gafun, gajac, gbfun, &
gbjac, guess
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: half = 0.5_nag_wp
Real (Kind=nag_wp), Parameter :: three = 3.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: two = 2.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: xsplit = 30.0_nag_wp
Integer, Parameter, Public :: 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), Public, Save :: el, en, s
Contains
Subroutine ffun(x,y,neq,m,f,iuser,ruser)
! .. 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) :: 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,iuser,ruser)
! .. 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)
! .. 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,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)
ga(2) = ya(1,1)
ga(3) = ya(2,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,1)
gb(2) = yb(2,0) - one
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(2,1,1) = one
dgady(3,2,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,1) = one
dgbdy(2,2,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)
! .. Local Scalars ..
Real (Kind=nag_wp) :: ex, expmx
! .. Intrinsic Procedures ..
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 d02txfe_mod, Only: el, en, ffun, fjac, gafun, gajac, gbfun, gbjac, &
guess, m1, m2, mmax, neq, nin, nlbc, nleft, nout, &
nrbc, nright, one, s, two, xsplit
Use nag_library, Only: d02tlf, d02tvf, d02txf, d02tyf, d02tzf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: dx, el_init, ermx, s_init, xx
Integer :: i, iermx, ifail, ijermx, j, licomm, &
lrcomm, mxmesh, ncol, ncont, nmesh
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: mesh(:), rcomm(:)
Real (Kind=nag_wp) :: ruser(1), tol(neq), y(neq,0:mmax-1)
Integer, Allocatable :: icomm(:), ipmesh(:)
Integer :: iuser(2), m(neq)
! .. Intrinsic Procedures ..
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)
Allocate (mesh(mxmesh),ipmesh(mxmesh))
! Read problem (initial) parameters
Read (nin,*) en, el_init, s_init
! Initialize 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
! Workspace query to get size of rcomm and icomm
ifail = 0
Call d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,ruser,0, &
iuser,2,ifail)
lrcomm = iuser(1)
licomm = iuser(2)
Allocate (rcomm(lrcomm),icomm(licomm))
! 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,rcomm, &
lrcomm,icomm,licomm,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 d02tlf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rcomm,icomm,iuser, &
ruser,ifail)
If (ifail/=0) Then
Exit cont
End If
! Extract mesh
ifail = 0
Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rcomm,icomm, &
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,rcomm,icomm,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,rcomm,icomm,ifail)
Write (nout,99998) xx*el, y(1,0), y(2,0)
End Do
! Select mesh for continuation and update continuation parameters.
If (j<ncont) Then
el = two*el
s = 0.6_nag_wp*s
nmesh = (nmesh+1)/2
ifail = 0
Call d02txf(mxmesh,nmesh,mesh,ipmesh,rcomm,icomm,ifail)
End If
End Do cont
99999 Format (/,' Solution on original interval:',/,6X,'x',8X,'f',10X,'g')
99998 Format (1X,F8.2,2(1X,F10.4))
99997 Format (/,/,' Tolerance = ',E8.1,' L = ',F8.3,' S =',F7.4)
99996 Format (/,' Used a mesh of ',I4,' points',/,' Maximum error = ',E10.2, &
' in interval ',I4,' for component ',I4)
End Program d02txfe