! D02MVF Example Program Text
! Mark 30.0 Release. NAG Copyright 2024.
Module d02mvfe_mod
! D02MVF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: jac, resid
! .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp
Integer, Parameter, Public :: iset = 1, neq = 6, nin = 5, nout = 6
Integer, Parameter, Public :: nrw = 50 + 4*neq
Integer, Parameter, Public :: nwkjac = neq*(neq+1)
Integer, Parameter, Public :: sdysav = 8
Integer, Parameter, Public :: ldysav = neq
Contains
Subroutine resid(neq,t,y,ydot,r,ires)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: t
Integer, Intent (Inout) :: ires
Integer, Intent (In) :: neq
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: r(neq)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
! .. Executable Statements ..
If (ires==-1) Then
r(1) = -ydot(1)
r(2) = -ydot(2)
r(3) = -ydot(3)
r(4) = -ydot(4)
r(5:6) = 0.0E0_nag_wp
Else
r(1) = y(3) - y(6)*y(1) - ydot(1)
r(2) = y(4) - y(6)*y(2) - ydot(2)
r(3) = -y(5)*y(1) - ydot(3)
r(4) = -y(5)*y(2) - one - ydot(4)
r(5) = y(1)*y(3) + y(2)*y(4)
r(6) = y(1)**2 + y(2)**2 - one
End If
Return
End Subroutine resid
Subroutine jac(neq,t,y,ydot,h,d,p)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: d, h, t
Integer, Intent (In) :: neq
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: p(neq,neq)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
! .. Local Scalars ..
Real (Kind=nag_wp) :: hxd
! .. Executable Statements ..
hxd = h*d
p(1,1) = (one+hxd*y(6))
p(1,3) = -hxd
p(1,6) = hxd*y(1)
p(2,2) = (one+hxd*y(6))
p(2,4) = -hxd
p(2,6) = hxd*y(2)
p(3,1) = hxd*y(5)
p(3,3) = one
p(3,5) = hxd*y(1)
p(4,2) = hxd*y(5)
p(4,4) = one
p(4,5) = hxd*y(2)
p(5,1) = -hxd*y(3)
p(5,2) = -hxd*y(4)
p(5,3) = -hxd*y(1)
p(5,4) = -hxd*y(2)
p(6,1) = -two*hxd*y(1)
p(6,2) = -two*hxd*y(2)
Return
End Subroutine jac
End Module d02mvfe_mod
Program d02mvfe
! D02MVF Example Main Program
! .. Use Statements ..
Use d02mvfe_mod, Only: iset, jac, ldysav, neq, nin, nout, nrw, nwkjac, &
one, resid, sdysav
Use nag_library, Only: d02mvf, d02nby, d02ngf, d02nsf, nag_wp, x04abf
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: h0, hmax, hmin, t, tcrit, tout
Integer :: i, ifail, itask, itol, itrace, &
maxord, maxstp, mxhnil, outchn
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: atol(:), rtol(:), rwork(:), &
wkjac(:), y(:), ydot(:), ysav(:,:)
Real (Kind=nag_wp) :: con(3)
Integer :: inform(23)
Logical :: lderiv(2)
! .. Executable Statements ..
Write (nout,*) 'D02MVF Example Program Results'
! Skip heading in data file
Read (nin,*)
! Allocations based on number of differential equations (neq)
Allocate (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq),ydot(neq), &
ysav(ldysav,sdysav))
! Read algorithmic parameters
Read (nin,*) maxord, maxstp, mxhnil
Read (nin,*) hmin, hmax, h0
Read (nin,*) rtol(1), atol(1)
Read (nin,*) itrace, itol
Read (nin,*) t, tout
Read (nin,*) y(1:neq)
Read (nin,*) itask
! Set initial derivatives and other values
outchn = nout
Call x04abf(iset,outchn)
con(1:3) = 0.0_nag_wp
ydot(1) = y(3) - y(6)*y(1)
ydot(2) = y(4) - y(6)*y(2)
ydot(3) = -y(5)*y(1)
ydot(4) = -y(5)*y(2) - one
ydot(5) = -3.0_nag_wp*y(4)
ydot(6) = 0.0_nag_wp
tcrit = tout
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02mvf(neq,sdysav,maxord,con,tcrit,hmin,hmax,h0,maxstp,mxhnil, &
'AVERAGE-L2',rwork,ifail)
Write (nout,*)
Write (nout,99999) 'Pendulum problem with relative tolerance', rtol(1)
Write (nout,99999) ' and absolute tolerance', atol(1)
Write (nout,99998)(i,i=1,neq)
Write (nout,99997) t, y(1:neq)
ifail = 0
Call d02nsf(neq,neq,'ANALYTIC',nwkjac,rwork,ifail)
lderiv(1) = .True.
lderiv(2) = .True.
ifail = 0
Call d02ngf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,resid, &
ysav,sdysav,jac,wkjac,nwkjac,d02nby,lderiv,itask,itrace,ifail)
Write (nout,99997) t, y(1:neq)
99999 Format (1X,A40,' =',1P,E8.1)
99998 Format (/,1X,' t ',3X,6(' y',I1))
99997 Format (1X,F7.4,2X,6(F8.4))
End Program d02mvfe