! D02EJF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02ejfe_mod ! Data for D02EJF example program ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. 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 :: n = 3, nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: h, xend Integer, Save :: k Contains Subroutine fcn(x,y,f) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: x ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: 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 (Out) :: 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 :: real ! .. Executable Statements .. Write (nout,99999) xsol, (y(j),j=1,n) xsol = xend - real(k,kind=nag_wp)*h k = k - 1 Return 99999 Format (1X,F8.2,3F13.5) End Subroutine output End Module d02ejfe_mod Program d02ejfe ! D02EJF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02ejf, d02ejw, d02ejx, d02ejy, nag_wp Use d02ejfe_mod, Only: fcn, g, h, k, n, nin, nout, output, pederv, xend ! .. 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)' k = kinit h = (xend-x)/real(k+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) Write (nout,*) ' Range too short for TOL' 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