dvindy Subroutine

private subroutine dvindy(me, t, k, yh, ldyh, dky, iflag)

dvindy computes interpolated values of the k-th derivative of the dependent variable vector y, and stores it in dky. this routine is called within the package with k = 0 and t = tout, but may also be called by the user for any k up to the current order. (see detailed instructions in the usage documentation.)

 the computed values in dky are gotten by interpolation using the
 nordsieck history array yh.  this array corresponds uniquely to a
 vector-valued polynomial of degree nqcur or less, and dky is set
 to the k-th derivative of this polynomial at t.
 the formula for dky is:
              q
  dky(i)  =  sum  c(j,k) * (t - tn)**(j-k) * h**(-j) * yh(i,j+1)
             j=k
 where  c(j,k) = j*(j-1)*...*(j-k+1), q = nqcur, tn = tcur, h = hcur.
 the quantities  nq = nqcur, l = nq+1, n, tn, and h are
 communicated by common.  the above sum is done in reverse order.
 iflag is returned negative if either k or t is out of bounds.

 discussion above and comments in driver explain all variables.

Type Bound

dvode_t

Arguments

Type IntentOptional Attributes Name
class(dvode_t), intent(inout) :: me
real(kind=wp), intent(in) :: t

value of independent variable where answers are desired (normally the same as the t last returned by dvode). for valid results, t must lie between tcur - hu and tcur. (see optional output for tcur and hu.)
integer, intent(in) :: k

integer order of the derivative desired. k must satisfy 0 <= k <= nqcur, where nqcur is the current order (see optional output). the capability corresponding to k = 0, i.e. computing y(t), is already provided by dvode directly. since nqcur >= 1, the first derivative dy/dt is always available with dvindy.
real(kind=wp), intent(in) :: yh(ldyh,*)

the history array yh
integer, intent(in) :: ldyh

column length of yh, equal to the initial value of neq.
real(kind=wp), intent(out) :: dky(*)

a real array of length neq containing the computed value of the k-th derivative of y(t).
integer, intent(out) :: iflag

integer flag, returned as 0 if k and t were legal, -1 if k was illegal, and -2 if t was illegal. on an error return, a message is also written.

Calls

proc~~dvindy~~CallsGraph proc~dvindy dvode_module::dvode_t%dvindy proc~dscal dvode_blas_module::dscal proc~dvindy->proc~dscal proc~xerrwd dvode_module::dvode_t%xerrwd proc~dvindy->proc~xerrwd proc~ixsav dvode_module::dvode_t%ixsav proc~xerrwd->proc~ixsav

Called by

proc~~dvindy~~CalledByGraph proc~dvindy dvode_module::dvode_t%dvindy proc~dvode dvode_module::dvode_t%dvode proc~dvode->proc~dvindy

Source Code

      subroutine dvindy(me,t,k,yh,ldyh,dky,iflag)

      class(dvode_t),intent(inout) :: me
      real(wp),intent(in) :: t !! value of independent variable where answers are desired
                               !! (normally the same as the t last returned by [[dvode]]).
                               !! for valid results, t must lie between tcur - hu and tcur.
                               !! (see optional output for tcur and hu.)
      integer,intent(in) :: k !! integer order of the derivative desired.  k must satisfy
                              !! 0 <= k <= nqcur, where nqcur is the current order
                              !! (see optional output).  the capability corresponding
                              !! to k = 0, i.e. computing y(t), is already provided
                              !! by [[dvode]] directly.  since nqcur >= 1, the first
                              !! derivative dy/dt is always available with [[dvindy]].
      integer,intent(in) :: ldyh !! column length of yh, equal to the initial value of neq.
      real(wp),intent(in) :: yh(ldyh,*) !! the history array yh
      real(wp),intent(out) :: dky(*) !! a real array of length neq containing the computed value
                                     !! of the k-th derivative of y(t).
      integer,intent(out) :: iflag !! integer flag, returned as 0 if k and t were legal,
                                   !! -1 if k was illegal, and -2 if t was illegal.
                                   !! on an error return, a message is also written.

      real(wp) :: c , r , s , tfuzz , tn1 , tp
      integer :: i , ic , j , jb , jb2 , jj , jj1 , jp1
      character(len=80) :: msg

      real(wp),parameter :: hun = 100.0_wp

      iflag = 0
      if ( k<0 .or. k>me%dat%nq ) then
         msg = 'dvindy-- k (=i1) illegal      '
         call me%xerrwd(msg,30,51,1,1,k,0,0,zero,zero)
         iflag = -1
         return
      else
         tfuzz = hun*me%dat%uround*sign(abs(me%dat%tn)+abs(me%dat%hu),me%dat%hu)
         tp = me%dat%tn - me%dat%hu - tfuzz
         tn1 = me%dat%tn + tfuzz
         if ( (t-tp)*(t-tn1)>zero ) then
            msg = 'dvindy-- t (=r1) illegal      '
            call me%xerrwd(msg,30,52,1,0,0,0,1,t,zero)
            msg =  '      t not in interval tcur - hu (= r1) to tcur (=r2)      '
            call me%xerrwd(msg,60,52,1,0,0,0,2,tp,me%dat%tn)
            iflag = -2
            return
         else
            s = (t-me%dat%tn)/me%dat%h
            ic = 1
            if ( k/=0 ) then
               jj1 = me%dat%l - k
               do jj = jj1 , me%dat%nq
                  ic = ic*jj
               enddo
            endif
            c = real(ic,wp)
            do i = 1 , me%dat%n
               dky(i) = c*yh(i,me%dat%l)
            enddo
            if ( k/=me%dat%nq ) then
               jb2 = me%dat%nq - k
               do jb = 1 , jb2
                  j = me%dat%nq - jb
                  jp1 = j + 1
                  ic = 1
                  if ( k/=0 ) then
                     jj1 = jp1 - k
                     do jj = jj1 , j
                        ic = ic*jj
                     enddo
                  endif
                  c = real(ic,wp)
                  do i = 1 , me%dat%n
                     dky(i) = c*yh(i,jp1) + s*dky(i)
                  enddo
               enddo
               if ( k==0 ) return
            endif
         endif
      endif
      r = me%dat%h**(-k)
      call dscal(me%dat%n,r,dky,1)

   end subroutine dvindy