! D02NEF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module d02nefe_mod
! D02NEF 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 :: jac1, jac2, res1, res2
! .. Parameters ..
Real (Kind=nag_wp), Parameter :: alpha = 0.04_nag_wp
Real (Kind=nag_wp), Parameter :: beta = 1.0E4_nag_wp
Real (Kind=nag_wp), Parameter :: gamma = 3.0E7_nag_wp
Real (Kind=nag_wp), Parameter :: two = 2.0_nag_wp
Integer, Parameter, Public :: ml = 1, mu = 2, neq1 = 3, neq2 = 1, &
nin = 5, nout = 6
Contains
Subroutine myjac1(neq,ml,mu,t,y,ydot,pd,cj)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: cj, t
Integer, Intent (In) :: ml, mu, neq
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: pd(2*ml+mu+1,neq)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
! .. Local Scalars ..
Integer :: md, ms
! .. Executable Statements ..
! Main diagonal pdfull(i,i), i=1,neq
md = mu + ml + 1
pd(md,1) = -alpha - cj
pd(md,2) = -beta*y(3) - two*gamma*y(2) - cj
pd(md,3) = -cj
! 1 Subdiagonal pdfull(i+1:i), i=1,neq-1
ms = md + 1
pd(ms,1) = alpha
pd(ms,2) = two*gamma*y(2)
! First superdiagonal pdfull(i-1,i), i=2, neq
ms = md - 1
pd(ms,2) = beta*y(3)
pd(ms,3) = -beta*y(2)
! Second superdiagonal pdfull(i-2,i), i=3, neq
ms = md - 2
pd(ms,3) = beta*y(2)
Return
End Subroutine myjac1
Subroutine myjac2(neq,t,y,ydot,pd,cj)
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: cj, t
Integer, Intent (In) :: neq
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: pd(neq*neq)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
! .. Intrinsic Procedures ..
Intrinsic :: exp
! .. Executable Statements ..
pd(1) = -two*y(1) + 0.1E0_nag_wp*t*y(1)*exp(y(1))
Return
End Subroutine myjac2
Subroutine res1(neq,t,y,ydot,r,ires,iuser,ruser)
! .. 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 (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
Integer, Intent (Inout) :: iuser(*)
! .. Executable Statements ..
r(1) = -alpha*y(1) + beta*y(2)*y(3) - ydot(1)
r(2) = alpha*y(1) - beta*y(2)*y(3) - gamma*y(2)*y(2) - ydot(2)
r(3) = gamma*y(2)*y(2) - ydot(3)
Return
End Subroutine res1
Subroutine jac1(neq,t,y,ydot,pd,cj,iuser,ruser)
! .. Use Statements ..
Use nag_library, Only: d02nez
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: cj, t
Integer, Intent (In) :: neq
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: pd(*), ruser(*)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Integer :: ijac, ml, mu
! .. Executable Statements ..
ml = iuser(1)
mu = iuser(2)
ijac = iuser(3)
If (ijac==1) Then
Call myjac1(neq,ml,mu,t,y,ydot,pd,cj)
Else
Call d02nez(neq,t,y,ydot,pd,cj,iuser,ruser)
End If
Return
End Subroutine jac1
Subroutine res2(neq,t,y,ydot,r,ires,iuser,ruser)
! .. 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 (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
Integer, Intent (Inout) :: iuser(*)
! .. Intrinsic Procedures ..
Intrinsic :: exp
! .. Executable Statements ..
r(1) = 4.0_nag_wp - y(1)**2 + t*0.1E0_nag_wp*exp(y(1))
Return
End Subroutine res2
Subroutine jac2(neq,t,y,ydot,pd,cj,iuser,ruser)
! .. Use Statements ..
Use nag_library, Only: d02nez
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: cj, t
Integer, Intent (In) :: neq
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: pd(*), ruser(*)
Real (Kind=nag_wp), Intent (In) :: y(neq), ydot(neq)
Integer, Intent (Inout) :: iuser(*)
! .. Local Scalars ..
Integer :: ijac
! .. Executable Statements ..
ijac = iuser(1)
If (ijac==1) Then
Call myjac2(neq,t,y,ydot,pd,cj)
Else
Call d02nez(neq,t,y,ydot,pd,cj,iuser,ruser)
End If
Return
End Subroutine jac2
End Module d02nefe_mod
Program d02nefe
! D02NEF Example Main Program
! .. Use Statements ..
Use d02nefe_mod, Only: nout
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Executable Statements ..
Write (nout,*) 'D02NEF Example Program Results'
Call ex1
Call ex2
Contains
Subroutine ex1
! .. Use Statements ..
Use d02nefe_mod, Only: jac1, ml, mu, neq1, nin, res1
Use nag_library, Only: d02mcf, d02mwf, d02nef, d02npf
! .. Local Scalars ..
Real (Kind=nag_wp) :: h0, hmax, t, tout
Integer :: i, ifail, ijac, itask, itol, j, &
lcom, licom, maxord, neq
Character (8) :: jceval
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: atol(:), com(:), rtol(:), y(:), &
ydot(:)
Real (Kind=nag_wp) :: ruser(1)
Integer, Allocatable :: icom(:)
Integer :: iuser(3)
! .. Executable Statements ..
Write (nout,*)
Write (nout,*) 'D02NEF Example 1'
! Skip heading in data file
Read (nin,*)
Read (nin,*) maxord
neq = neq1
lcom = 40 + (maxord+4)*neq + (2*ml+mu+1)*neq + 2*(neq/(ml+mu+1)+1)
licom = 50 + neq
Allocate (atol(neq),com(lcom),rtol(neq),y(neq),ydot(neq),icom(licom))
Read (nin,*) ijac, itol
Read (nin,*) rtol(1:neq)
Read (nin,*) atol(1:neq)
Read (nin,*) ydot(1:neq)
If (ijac==1) Then
jceval = 'Analytic'
Else
jceval = 'Numeric'
End If
! Set initial values
Read (nin,*) y(1:neq)
! Initialize the problem, specifying that the Jacobian is to be
! evaluated analytically using the provided routine jac.
Read (nin,*) hmax, h0
Read (nin,*) t, tout
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02mwf(neq,maxord,jceval,hmax,h0,itol,icom,licom,com,lcom,ifail)
! Specify that the Jacobian is banded.
ifail = 0
Call d02npf(neq,ml,mu,icom,licom,ifail)
! Use the iuser array to pass the band dimensions through to jac.
! An alternative would be to hard code values for ml and mu in jac.
iuser(1) = ml
iuser(2) = mu
iuser(3) = ijac
Write (nout,99999)(i,i=1,neq)
Write (nout,99998) t, (y(i),i=1,neq)
itask = 0
! Obtain the solution at 5 equally spaced values of T.
loop: Do j = 1, 5
ifail = -1
Call d02nef(neq,t,tout,y,ydot,rtol,atol,itask,res1,jac1,icom,com, &
lcom,iuser,ruser,ifail)
Write (nout,99998) t, (y(i),i=1,neq)
If (ifail/=0) Then
Write (nout,99997) ifail
Exit loop
End If
tout = tout + 0.02_nag_wp
Call d02mcf(icom)
End Do loop
Write (nout,*)
Write (nout,99996) itask
99999 Format (/,1X,' t ',5X,3(' Y(',I1,') '))
99998 Format (1X,F8.4,3X,3(F12.6))
99997 Format (1X,' ** D02NEF returned with IFAIL = ',I5)
99996 Format (1X,'The integrator completed task, ITASK = ',I4)
End Subroutine ex1
Subroutine ex2
! .. Use Statements ..
Use d02nefe_mod, Only: jac2, neq2, nin, res2
Use nag_library, Only: d02mcf, d02mwf, d02nef
! .. Local Scalars ..
Real (Kind=nag_wp) :: h0, hmax, t, tout
Integer :: i, ifail, ijac, itask, itol, j, &
lcom, licom, maxord, neq
Character (8) :: jceval
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: atol(:), com(:), rtol(:), y(:), &
ydot(:)
Real (Kind=nag_wp) :: ruser(1)
Integer, Allocatable :: icom(:)
Integer :: iuser(1)
! .. Executable Statements ..
Write (nout,*)
Write (nout,*) 'D02NEF Example 2'
Write (nout,*)
Read (nin,*)
Read (nin,*) maxord
neq = neq2
lcom = 40 + (maxord+4)*neq + neq*neq
licom = 50 + neq
Allocate (atol(neq),com(lcom),rtol(neq),y(neq),ydot(neq),icom(licom))
Read (nin,*) ijac, itol
Read (nin,*) rtol(1:neq)
Read (nin,*) atol(1:neq)
Read (nin,*) ydot(1:neq)
If (ijac==1) Then
jceval = 'Analytic'
Else
jceval = 'Numeric'
End If
! Initialize the problem, specifying that the Jacobian is to be
! evaluated analytically using the provided routine jac.
Read (nin,*) y(1:neq)
Read (nin,*) hmax, h0
Read (nin,*) t, tout
ifail = 0
Call d02mwf(neq,maxord,jceval,hmax,h0,itol,icom,licom,com,lcom,ifail)
! Use the iuser array to pass whether numerical or analytic Jacobian
! is to be used.
iuser(1) = ijac
Write (nout,99999)(i,i=1,neq)
Write (nout,99998) t, y(1:neq)
itask = 0
! Obtain the solution at 5 equally spaced values of t.
loop: Do j = 1, 5
ifail = -1
Call d02nef(neq,t,tout,y,ydot,rtol,atol,itask,res2,jac2,icom,com, &
lcom,iuser,ruser,ifail)
Write (nout,99998) t, y(1:neq)
If (ifail/=0) Then
Write (nout,99997) ifail
Exit loop
End If
tout = tout + 0.2_nag_wp
Call d02mcf(icom)
End Do loop
Write (nout,*)
Write (nout,99996) itask
99999 Format (/,1X,' t y(',I1,')')
99998 Format (1X,F8.4,3X,3(F12.6))
99997 Format (1X,' ** D02NEF returned with IFAIL = ',I5)
99996 Format (1X,'The integrator completed task, ITASK = ',I4)
End Subroutine ex2
End Program d02nefe