Here’s an example of using dvindy for dense output.
Example
program droplet
implicit none
integer, parameter :: dp = kind(1.0d0)
real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp)
integer, parameter :: neq = 3
integer :: itol, itask, istate, iopt, mf, iflag, n, ngrid
real(dp) :: rtol, s, send, ds
real(dp), dimension(neq) :: atol, y
real(dp), allocatable :: ygrid(:,:), sgrid(:)
real(dp), allocatable :: rwork(:)
integer, allocatable :: iwork(:)
integer :: lrw, liw
real(dp) :: rpar(2)
integer :: idum(1)
external :: dvode, dvindy
rpar(1) = 3.3888_dp
rpar(2) = 0.8414_dp
s = 0.0_dp
send = 2.0_dp
ngrid = 41
allocate(sgrid(ngrid),ygrid(neq,ngrid))
sgrid(1) = 0.0_dp
ygrid(:,1) = [0.0_dp, 0.0_dp, 0.0_dp]
ds = (send - s)/real(ngrid - 1,dp)
do n = 2, ngrid
sgrid(n) = (n - 1)*ds
end do
n = 2
ds = 0.00001_dp
s = ds
y(1) = ds
y(2) = cos(y(1))*ds
y(3) = sin(y(1))*ds
itol = 2
rtol = 1.0e-7_dp
atol(1) = 1.0e-6_dp
atol(2) = 1.0e-6_dp
atol(3) = 1.0e-6_dp
mf = 22
lrw = len_rwork(neq,mf)
liw = len_iwork(neq,mf)
allocate(rwork(lrw), iwork(liw))
istate = 1
itask = 2
iopt = 0
integrate: do
call dvode(young_laplace, neq, y, s, send, &
itol, rtol, atol, itask, istate, iopt, &
rwork, lrw, iwork, liw, &
young_laplace_jacobian, &
mf, rpar, idum)
! Automatic stepping
write(*,'(4(ES12.4))') s, y
if (istate < 0) then
write(*,*) ' Failed at s = ',s,' with istate = ',istate
error stop
end if
! Interpolation
interpolate: do while(s >= sgrid(n))
call dvindy (sgrid(n), 0, rwork(21), neq, ygrid(:,n), iflag)
n = n + 1
if (n > ngrid) exit interpolate
end do interpolate
if (s >= send) then
exit integrate
end if
end do integrate
write(*,*)
write(*,*)
do n = 1, ngrid
write(*,'(4(ES12.4))') sgrid(n), ygrid(:,n)
end do
contains
subroutine young_laplace(neq, t, y, ydot, rpar, ipar)
integer, intent(in) :: neq
real(dp), intent(in) :: t, y(neq)
real(dp), intent(out) :: ydot(neq)
real(dp), intent(in) :: rpar(*)
integer, intent(in) :: ipar(*)
associate(pL => rpar(1), drho => rpar(2))
ydot(1) = pL - drho*y(3) - sin(y(1))/y(2)
ydot(2) = cos(y(1))
ydot(3) = sin(y(1))
end associate
end subroutine
subroutine young_laplace_jacobian(neq, t, y, ml, mu, pd, nrpd, rpar, ipar)
integer, intent(in) :: neq, nrpd, ml, mu
real(dp), intent(out) :: pd(nrpd,neq)
real(dp), intent(in) :: t, y(neq), rpar(*)
integer, intent(in) :: ipar(*)
! dummy routine
associate(pL => rpar(1), drho => rpar(2))
pd(1,1) = -cos(y(1))/y(2)
pd(2,1) = -sin(y(1))
pd(3,1) = cos(y(1))
pd(1,2) = sin(y(1))/(y(2))**2
pd(2,2) = 0.0_dp
pd(3,2) = 0.0_dp
pd(1,3) = -drho
pd(2,3) = 0.0_dp
pd(3,3) = 0.0_dp
end associate
end subroutine
function len_rwork(neq, mf, mlu) result(lrw)
integer, intent(in) :: neq, mf
integer, intent(in), optional :: mlu(2)
! Upper and lower bandwidth in case of banded Jacobian
integer :: lrw
! RWORK = Real work array of length at least:
! 20 + 16*NEQ for MF = 10,
! 22 + 9*NEQ + 2*NEQ**2 for MF = 21 or 22,
! 22 + 11*NEQ + (3*ML + 2*MU)*NEQ for MF = 24 or 25.
select case(mf)
case(10)
lrw = 20 + 17*neq
case(21,22)
lrw = 22 + 9*neq + 2*neq**2
case(24,25)
if (present(mlu)) then
lrw = 22 + 11*neq + (3*mlu(1) + 2*mlu(2))*neq
else
write(*,*) "len_rwork: Missing array mlu containing bandwidths [ml, mu]"
error stop
end if
case default
! Wrong value for mf
lrw = -1
end select
end function
function len_iwork(neq, mf) result(liw)
integer, intent(in) :: neq, mf
integer :: liw
! IWORK = Integer work array of length at least:
! 30 for MF = 10,
! 30 + NEQ for MF = 21, 22, 24, or 25.
! If MF = 24 or 25, input in IWORK(1),IWORK(2) the lower
! and upper half-bandwidths ML,MU.
select case (mf)
case (10)
liw = 30
case (21, 22, 24, 25)
liw = 30 + neq
case default
! Wrong value for mf
liw = -1
end select
end function
end program
The system solved are the Young-Laplace equations for a pendant drop in a cylindrical coordinate system:
An equivalent system of equations (but with a slightly different set of parameters) was first solved in 1883 by Francis Bashforth and John Couch Adams in An attempt to test the theories of capillary action by comparing the theoretical and measured forms of drops of fluid - GDZ. Notably, for this occasion Adams derived the family of explicit multi-step methods, now called Adams-Bashforth methods.
In VODE, when using mf = 10 the solver uses an implicit method from Adams. Here is plot of the interpolated (r,z) pairs super-imposed on the automatic step control:
The resulting profile traces half of a pendant drop. The interpolated points are spaced equidistantly in terms of the arc length s.
By fitting the theoretical drop shape to an experimental shape profile, the surface tension of a liquid can be determined. The method is called pendant drop tensiometry.