! D02LAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 ! .. 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 Procedures .. 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