! D02CJF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d02cjfe_mod ! Data for D02CJF example program ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: n = 3, nin = 5, nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: h, xend Integer, Save :: k ! n: number of differential equations Contains 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 Subroutine fcn(x,y,f) ! .. Parameters .. Real (Kind=nag_wp), Parameter :: alpha = -0.032E0_nag_wp Real (Kind=nag_wp), Parameter :: beta = -0.02E0_nag_wp ! .. 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(*) ! .. Intrinsic Procedures .. Intrinsic :: cos, tan ! .. Executable Statements .. f(1) = tan(y(3)) f(2) = alpha*tan(y(3))/y(2) + beta*y(2)/cos(y(3)) f(3) = alpha/y(2)**2 Return End Subroutine fcn 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) Return End Function g End Module d02cjfe_mod Program d02cjfe ! D02CJF Example Main Program ! .. Use Statements .. Use nag_library, Only: d02cjf, d02cjw, d02cjx, nag_wp Use d02cjfe_mod, Only: fcn, g, h, k, n, nin, nout, output, 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,*) 'D02CJF Example Program Results' iw = 21*n + 28 Allocate (w(iw),y(n),yinit(n)) ! Skip heading in data file Read (nin,*) ! xinit: initial x value, xend: final x value. Read (nin,*) xinit Read (nin,*) xend Read (nin,*) yinit(1:n) Read (nin,*) kinit Do icase = 1, 4 Write (nout,*) Select Case (icase) Case (1) Write (nout,99995) icase, 'intermediate output, root-finding' Case (2) Write (nout,99995) icase, 'no intermediate output, root-finding' Case (3) Write (nout,99995) icase, 'intermediate output, no root-finding' Case (4) Write (nout,99995) icase, & 'no intermediate output, no root-finding ( integrate to XEND)' End Select Do j = 4, 5 tol = 10.0E0_nag_wp**(-j) Write (nout,*) Write (nout,99999) ' Calculation with TOL =', tol x = xinit y(1:n) = yinit(1:n) If (icase/=2) Then Write (nout,*) ' X Y(1) Y(2) Y(3)' k = kinit h = (xend-x)/real(k+1,kind=nag_wp) End If ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Select Case (icase) Case (1) Call d02cjf(x,xend,n,y,fcn,tol,'Default',output,g,w,ifail) Write (nout,99998) ' Root of Y(1) = 0.0 at', x Write (nout,99997) ' Solution is', (y(i),i=1,n) Case (2) Call d02cjf(x,xend,n,y,fcn,tol,'Default',d02cjx,g,w,ifail) Write (nout,99998) ' Root of Y(1) = 0.0 at', x Write (nout,99997) ' Solution is', (y(i),i=1,n) Case (3) Call d02cjf(x,xend,n,y,fcn,tol,'Default',output,d02cjw,w,ifail) Case (4) Write (nout,99996) x, (y(i),i=1,n) Call d02cjf(x,xend,n,y,fcn,tol,'Default',d02cjx,d02cjw,w,ifail) Write (nout,99996) x, (y(i),i=1,n) End Select End Do If (icase<4) 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,': ',A) End Program d02cjfe