! D02HBF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02hbfe_mod ! Data for D02HBF example programs ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: iset = 1, nin = 5, nout = 6 CONTAINS SUBROUTINE fcn1(x,y,f,p) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: f(*) REAL (KIND=nag_wp), INTENT (IN) :: p(*), y(*) ! .. Executable Statements .. f(1) = y(2) f(2) = (y(1)**3-y(2))/(2.0E0_nag_wp*x) RETURN END SUBROUTINE fcn1 SUBROUTINE range1(a,b,p) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: a, b ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: p(*) ! .. Executable Statements .. a = 0.1E0_nag_wp b = 16.0E0_nag_wp RETURN END SUBROUTINE range1 SUBROUTINE bc1(g1,g2,p) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: g1(*), g2(*) REAL (KIND=nag_wp), INTENT (IN) :: p(*) ! .. Local Scalars .. REAL (KIND=nag_wp) :: z ! .. Intrinsic Functions .. INTRINSIC sqrt ! .. Executable Statements .. z = 0.1E0_nag_wp g1(1) = 0.1E0_nag_wp + p(1)*sqrt(z)*0.1E0_nag_wp + 0.01E0_nag_wp*z g1(2) = p(1)*0.05E0_nag_wp/sqrt(z) + 0.01E0_nag_wp g2(1) = 1.0E0_nag_wp/6.0E0_nag_wp g2(2) = p(2) RETURN END SUBROUTINE bc1 SUBROUTINE fcn2(x,y,f,p) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: f(*) REAL (KIND=nag_wp), INTENT (IN) :: p(*), y(*) ! .. Intrinsic Functions .. INTRINSIC cos, tan ! .. Executable Statements .. f(1) = tan(y(3)) f(2) = -p(1)*tan(y(3))/y(2) - 0.00002E0_nag_wp*y(2)/cos(y(3)) f(3) = -p(1)/y(2)**2 RETURN END SUBROUTINE fcn2 SUBROUTINE range2(a,b,p) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: a, b ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: p(*) ! .. Executable Statements .. a = 0.0E0_nag_wp b = p(2) RETURN END SUBROUTINE range2 SUBROUTINE bc2(g1,g2,p) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: g1(*), g2(*) REAL (KIND=nag_wp), INTENT (IN) :: p(*) ! .. Executable Statements .. g1(1) = 0.0E0_nag_wp g1(2) = 500.0E0_nag_wp g1(3) = 0.5E0_nag_wp g2(1) = 0.0E0_nag_wp g2(2) = 450.0E0_nag_wp g2(3) = p(3) RETURN END SUBROUTINE bc2 END MODULE d02hbfe_mod PROGRAM d02hbfe ! D02HBF Example Main Program ! .. Use Statements .. USE d02hbfe_mod, ONLY : nout ! .. Implicit None Statement .. IMPLICIT NONE ! .. Executable Statements .. WRITE (nout,*) 'D02HBF Example Program Results' CALL ex1 CALL ex2 CONTAINS SUBROUTINE ex1 ! .. Use Statements .. USE nag_library, ONLY : d02hbf, nag_wp, x04abf USE d02hbfe_mod, ONLY : bc1, fcn1, iset, nin, range1 ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: h, x, x1, xh INTEGER :: i, ifail, m1, n, n1, outchn, & sdw ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: e(:), p(:), pe(:), soln(:,:), & w(:,:) ! .. Intrinsic Functions .. INTRINSIC real ! .. Executable Statements .. ! Skip heading in data file READ (nin,*) ! m1: controls exit values, n: number of differential equations, ! n1: number of parameters. READ (nin,*) m1, n, n1 sdw = 3*n + 14 + 11 ALLOCATE (e(n),p(n1),pe(n1),soln(n,m1),w(n,sdw)) WRITE (nout,*) outchn = nout WRITE (nout,*) CALL x04abf(iset,outchn) ! p: estimates for the parameters p, e: bound on the local error. READ (nin,*) p(1:n1) READ (nin,*) pe(1:n1) READ (nin,*) e(1:n) WRITE (nout,*) 'Case 1' WRITE (nout,*) ! ifail: behaviour on error exit ! =1 for quiet-soft exit ! * Set ifail to 111 to obtain monitoring information * ifail = 1 CALL d02hbf(p,n1,pe,e,n,soln,m1,fcn1,bc1,range1,w,sdw,ifail) IF (ifail==0) THEN WRITE (nout,*) 'Final parameters' WRITE (nout,99999) (p(i),i=1,n1) WRITE (nout,*) WRITE (nout,*) 'Final solution' WRITE (nout,*) 'X-value Components of solution' CALL range1(x,x1,p) h = (x1-x)/real(m1-1,kind=nag_wp) xh = x DO i = 1, m1 WRITE (nout,99998) xh, soln(1:n,i) xh = xh + h END DO ELSE WRITE (nout,99996) ifail IF (ifail>1 .AND. ifail<=5) THEN WRITE (nout,99997) w(1,2), (w(i,1),i=1,n) END IF END IF RETURN 99999 FORMAT (1X,1P,3E15.3) 99998 FORMAT (1X,F7.2,2F13.4) 99997 FORMAT (/1X,'W(1,2) = ',F9.4,' W(.,1) = ',10E10.3) 99996 FORMAT (1X/1X,' ** D02HBF returned with IFAIL = ',I5) END SUBROUTINE ex1 SUBROUTINE ex2 ! .. Use Statements .. USE nag_library, ONLY : d02hbf, nag_wp, x04abf USE d02hbfe_mod, ONLY : bc2, fcn2, iset, nin, range2 ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: h, x, x1, xh INTEGER :: i, ifail, m1, n, n1, outchn, & sdw ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: e(:), p(:), pe(:), soln(:,:), & w(:,:) ! .. Intrinsic Functions .. INTRINSIC real ! .. Executable Statements .. READ (nin,*) ! m1: controls exit values, n: number of differential equations, ! n1: number of parameters. READ (nin,*) m1, n, n1 sdw = 3*n + 14 + 11 ALLOCATE (e(n),p(n1),pe(n1),soln(n,m1),w(n,sdw)) outchn = nout CALL x04abf(iset,outchn) ! p: estimates for the parameters p, e: bound on the local error. READ (nin,*) p(1:n1) READ (nin,*) pe(1:n1) READ (nin,*) e(1:n) WRITE (nout,*) WRITE (nout,*) WRITE (nout,*) 'Case 2' WRITE (nout,*) ! ifail: behaviour on error exit ! =1 for quiet-soft exit ! * Set ifail to 111 to obtain monitoring information * ifail = 1 CALL d02hbf(p,n1,pe,e,n,soln,m1,fcn2,bc2,range2,w,sdw,ifail) IF (ifail==0) THEN WRITE (nout,*) 'Final parameters' WRITE (nout,99999) (p(i),i=1,n1) WRITE (nout,*) WRITE (nout,*) 'Final solution' WRITE (nout,*) 'X-value Components of solution' CALL range2(x,x1,p) h = (x1-x)/real(m1-1,kind=nag_wp) xh = x DO i = 1, m1 WRITE (nout,99998) xh, soln(1:n,i) xh = xh + h END DO ELSE WRITE (nout,99996) ifail IF (ifail>1 .AND. ifail<=5) THEN WRITE (nout,99997) w(1,2), (w(i,1),i=1,n) END IF END IF RETURN 99999 FORMAT (1X,1P,3E15.3) 99998 FORMAT (1X,F7.0,2F13.1,F13.3) 99997 FORMAT (/1X,'W(1,2) = ',F9.4,' W(.,1) = ',10E10.3) 99996 FORMAT (1X/1X,' ** D02HBF returned with IFAIL = ',I5) END SUBROUTINE ex2 END PROGRAM d02hbfe