! D02LAF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d02lafe_mod ! D02LAF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. REAL (KIND=nag_wp), PARAMETER :: zero = 0.0_nag_wp INTEGER, PARAMETER :: neq = 2, nin = 5, nout = 6 INTEGER, PARAMETER :: lrwork = 16 + 20*neq CONTAINS SUBROUTINE fcn(neq,t,y,f) ! Derivatives for two body problem in y'' = f(t,y) form ! .. Implicit None Statement .. IMPLICIT NONE ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: t INTEGER, INTENT (IN) :: neq ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (OUT) :: f(neq) REAL (KIND=nag_wp), INTENT (IN) :: y(neq) ! .. Local Scalars .. REAL (KIND=nag_wp) :: r ! .. Intrinsic Functions .. INTRINSIC sqrt ! .. Executable Statements .. r = sqrt(y(1)**2+y(2)**2)**3 f(1) = -y(1)/r f(2) = -y(2)/r RETURN END SUBROUTINE fcn END MODULE d02lafe_mod PROGRAM d02lafe ! D02LAF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : d02laf, d02lxf, d02lyf, d02lzf, nag_wp USE d02lafe_mod, ONLY : fcn, lrwork, neq, nin, nout, zero ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: h, hnext, hstart, hused, t, & tend, tinc, tnext, tol, tstart INTEGER :: i, ifail, itol, maxstp, natt, & nfail, nsucc, nwant LOGICAL :: high, onestp, start ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: rwork(:), thres(:), thresp(:), & y(:), ydp(:), yinit(:), yp(:), & ypinit(:), ypwant(:), ywant(:) ! .. Executable Statements .. WRITE (nout,*) 'D02LAF Example Program Results' ! Skip heading in data file READ (nin,*) ! neq: number of second-order ordinary differential equations READ (nin,*) nwant ALLOCATE (rwork(lrwork),thres(neq),thresp(neq),y(neq),ydp(neq), & yinit(neq),yp(neq),ypinit(neq),ypwant(nwant),ywant(nwant)) READ (nin,*) high, onestp READ (nin,*) tinc ! Initial conditions READ (nin,*) tstart, tend READ (nin,*) yinit(1:neq) READ (nin,*) ypinit(1:neq) LOOP1: DO itol = 4, 5 tol = 10.0_nag_wp**(-itol) WRITE (nout,*) ! Call D02LXF with default THRES,THRESP,MAXSTP and H thres(1) = zero thresp(1) = zero h = zero maxstp = 0 start = .TRUE. ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 CALL d02lxf(neq,h,tol,thres,thresp,maxstp,start,onestp,high,rwork, & lrwork,ifail) WRITE (nout,99999) 'Calculation with TOL = ', tol WRITE (nout,99995) (i,i=1,neq) ! Set initial values y(1:neq) = yinit(1:neq) yp(1:neq) = ypinit(1:neq) t = tstart tnext = t + tinc WRITE (nout,99998) t, y(1:neq) ! Loop point for onestep mode LOOP2: DO ifail = -1 CALL d02laf(fcn,neq,t,tend,y,yp,ydp,rwork,lrwork,ifail) IF (ifail>0) THEN WRITE (nout,99997) ifail, t EXIT LOOP1 END IF ! Loop point for interpolation DO WHILE (tnext<=t) ifail = 0 CALL d02lzf(neq,t,y,yp,neq,tnext,ywant,ypwant,rwork,lrwork, & ifail) WRITE (nout,99998) tnext, ywant(1:neq) tnext = tnext + tinc END DO IF (t>=tend) EXIT LOOP2 END DO LOOP2 ifail = 0 CALL d02lyf(neq,hnext,hused,hstart,nsucc,nfail,natt,thres,thresp, & rwork,lrwork,ifail) WRITE (nout,*) WRITE (nout,99996) ' Number of successful steps = ', nsucc WRITE (nout,99996) ' Number of failed steps = ', nfail END DO LOOP1 99999 FORMAT (1X,A,1P,E9.1) 99998 FORMAT (1X,F5.1,2(2X,F9.5)) 99997 FORMAT (/1X,'D02LAF returned with IFAIL = ',I2,' at T = ',1P,E10.3) 99996 FORMAT (1X,A,I5) 99995 FORMAT (/' T ',2(' Y(',I1,') ')) END PROGRAM d02lafe