! D02EJF Example Program Text
! Mark 26.1 Release. NAG Copyright 2017.
Module d02ejfe_mod
! Data for D02EJF example program
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: fcn, g, output, pederv
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: alpha = 0.04_nag_wp
Real (Kind=nag_wp), Parameter :: beta = 1.0E4_nag_wp
Real (Kind=nag_wp), Parameter :: gamma = 3.0E7_nag_wp
Real (Kind=nag_wp), Parameter :: zero = 0.0_nag_wp
Integer, Parameter, Public :: n = 3, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: h, xend
Contains
Subroutine fcn(x,y,f)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: f(*)
Real (Kind=nag_wp), Intent (In) :: y(*)
! .. Executable Statements ..
f(1) = -alpha*y(1) + beta*y(2)*y(3)
f(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2)
f(3) = gamma*y(2)*y(2)
Return
End Subroutine fcn
Subroutine pederv(x,y,pw)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: pw(*)
Real (Kind=nag_wp), Intent (In) :: y(*)
! .. Executable Statements ..
pw(1) = -alpha
pw(2) = alpha
pw(3) = zero
pw(4) = beta*y(3)
pw(5) = -beta*y(3) - 2.0_nag_wp*gamma*y(2)
pw(6) = 2.0_nag_wp*gamma*y(2)
pw(7) = beta*y(2)
pw(8) = -beta*y(2)
pw(9) = zero
Return
End Subroutine pederv
Function g(x,y)
! .. Function Return Value ..
Real (Kind=nag_wp) :: g
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: y(*)
! .. Executable Statements ..
g = y(1) - 0.9E0_nag_wp
Return
End Function g
Subroutine output(xsol,y)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: xsol
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: y(*)
! .. Local Scalars ..
Integer :: j
! .. Intrinsic Procedures ..
Intrinsic :: abs
! .. Executable Statements ..
Write (nout,99999) xsol, (y(j),j=1,n)
xsol = xsol + h
! Make sure we exactly hit xsol = xend
If (abs(xsol-xend)<h/4.0E0_nag_wp) Then
xsol = xend
End If
Return
99999 Format (1X,F8.2,3F13.5)
End Subroutine output
End Module d02ejfe_mod
Program d02ejfe
! D02EJF Example Main Program
! .. Use Statements ..
Use d02ejfe_mod, Only: fcn, g, h, n, nin, nout, output, pederv, xend
Use nag_library, Only: d02ejf, d02ejw, d02ejx, d02ejy, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: tol, x, xinit
Integer :: i, icase, ifail, iw, j, kinit
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: w(:), y(:), yinit(:)
! .. Intrinsic Procedures ..
Intrinsic :: real
! .. Executable Statements ..
Write (nout,*) 'D02EJF Example Program Results'
iw = (12+n)*n + 50
Allocate (w(iw),y(n),yinit(n))
! Skip heading in data file
Read (nin,*)
! xinit: initial x value, xend: final x value
! y: initial solution values
Read (nin,*) xinit, xend
Read (nin,*) yinit(1:n)
Read (nin,*) kinit
Do icase = 1, 5
If (icase/=2) Then
Write (nout,99995) icase, 'Jacobian internally'
Else
Write (nout,99995) icase, 'Jacobian by PEDERV'
End If
Select Case (icase)
Case (1,2)
Write (nout,99994) 'intermediate output, root-finding'
Case (3)
Write (nout,99994) 'no intermediate output, root-finding'
Case (4)
Write (nout,99994) 'intermediate output, no root-finding'
Case (5)
Write (nout,99994) &
'no intermediate output, no root-finding (integrate to XEND)'
End Select
Do j = 3, 4
tol = 10.0E0_nag_wp**(-j)
Write (nout,99999) ' Calculation with TOL =', tol
x = xinit
y(1:n) = yinit(1:n)
If (icase/=3) Then
Write (nout,*) ' X Y(1) Y(2) Y(3)'
h = (xend-x)/real(kinit+1,kind=nag_wp)
End If
ifail = 0
Select Case (icase)
Case (1)
Call d02ejf(x,xend,n,y,fcn,d02ejy,tol,'Default',output,g,w,iw, &
ifail)
Write (nout,99998) ' Root of Y(1)-0.9 at', x
Write (nout,99997) ' Solution is', (y(i),i=1,n)
Case (2)
Call d02ejf(x,xend,n,y,fcn,pederv,tol,'Default',output,g,w,iw, &
ifail)
Write (nout,99998) ' Root of Y(1)-0.9 at', x
Write (nout,99997) ' Solution is', (y(i),i=1,n)
Case (3)
Call d02ejf(x,xend,n,y,fcn,d02ejy,tol,'Default',d02ejx,g,w,iw, &
ifail)
Write (nout,99998) ' Root of Y(1)-0.9 at', x
Write (nout,99997) ' Solution is', (y(i),i=1,n)
Case (4)
ifail = 0
Call d02ejf(x,xend,n,y,fcn,d02ejy,tol,'Default',output,d02ejw,w, &
iw,ifail)
Case (5)
Write (nout,99996) x, (y(i),i=1,n)
Call d02ejf(x,xend,n,y,fcn,d02ejy,tol,'Default',d02ejx,d02ejw,w, &
iw,ifail)
Write (nout,99996) x, (y(i),i=1,n)
End Select
If (tol<0.0E0_nag_wp) Then
Write (nout,*) ' Range too short for TOL'
End If
End Do
If (icase<5) Then
Write (nout,*)
End If
End Do
99999 Format (/,1X,A,E8.1)
99998 Format (1X,A,F7.3)
99997 Format (1X,A,3F13.5)
99996 Format (1X,F8.2,3F13.5)
99995 Format (/,1X,'Case ',I1,': calculating ',A,',')
99994 Format (8X,A)
End Program d02ejfe