dintp Subroutine

private subroutine dintp(x, y, xout, yout, ypout, neqn, kold, phi, ivc, iv, kgi, gi, alpha, og, ow, ox, oy)

approximate the solution at xout by evaluating the polynomial computed in dsteps at xout. must be used in conjunction with dsteps.

the methods in subroutine dsteps approximate the solution near x by a polynomial. subroutine dintp approximates the solution at xout by evaluating the polynomial there. information defining this polynomial is passed from dsteps so dintp cannot be used alone.

References

  • L. F. Shampine, M. K. Gordon, "Computer solution of ordinary differential equations, the initial value problem", W. H. Freeman and Company, 1975.
  • H. A. Watts, "A smoother interpolant for DE/STEP, INTRP and DEABM: II", Report SAND84-0293, Sandia Laboratories, 1984.

History

  • 840201 date written -- watts, h. a., (snla)
  • 890831 modified array declarations. (wrb)
  • 890831 revision date from version 3.2
  • 891214 prologue converted to version 4.0 format. (bab)
  • 920501 reformatted the references section. (wrb)
  • December, 2015 : Refactored this routine (jw)

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: x
real(kind=wp), intent(in), dimension(neqn) :: y
real(kind=wp), intent(in) :: xout

point at which solution is desired
real(kind=wp), intent(out), dimension(neqn) :: yout

solution at xout
real(kind=wp), intent(out), dimension(neqn) :: ypout

derivative of solution at xout
integer, intent(in) :: neqn
integer, intent(in) :: kold
real(kind=wp), intent(in), dimension(neqn, 16) :: phi
integer, intent(in) :: ivc
integer, intent(in), dimension(10) :: iv
integer, intent(in) :: kgi
real(kind=wp), intent(in), dimension(11) :: gi
real(kind=wp), intent(in), dimension(12) :: alpha
real(kind=wp), intent(in), dimension(13) :: og
real(kind=wp), intent(in), dimension(12) :: ow
real(kind=wp), intent(in) :: ox
real(kind=wp), intent(in), dimension(neqn) :: oy

Called by

proc~~dintp~~CalledByGraph proc~dintp dintp proc~ddeabm_interp ddeabm_class%ddeabm_interp proc~ddeabm_interp->proc~dintp proc~ddes ddeabm_class%ddes proc~ddes->proc~dintp proc~ddeabm ddeabm_class%ddeabm proc~ddeabm->proc~ddes proc~ddeabm_with_event_wrapper ddeabm_with_event_class%ddeabm_with_event_wrapper proc~ddeabm_with_event_wrapper->proc~ddeabm_interp proc~ddeabm_with_event_wrapper->proc~ddeabm proc~ddeabm_with_event_wrapper_vec ddeabm_with_event_class_vec%ddeabm_with_event_wrapper_vec proc~ddeabm_with_event_wrapper_vec->proc~ddeabm proc~ddeabm_wrapper ddeabm_class%ddeabm_wrapper proc~ddeabm_wrapper->proc~ddeabm

Source Code

   subroutine dintp(x, y, xout, yout, ypout, neqn, kold, phi, ivc, iv, kgi, gi, alpha, og, ow, ox, oy)

      implicit none

      integer, intent(in)                     :: neqn
      real(wp), intent(in)                    :: x
      real(wp), dimension(neqn), intent(in)    :: y
      real(wp), intent(in)                    :: xout   !! point at which solution is desired
      real(wp), dimension(neqn), intent(out)   :: yout   !! solution at xout
      real(wp), dimension(neqn), intent(out)   :: ypout  !! derivative of solution at xout
      integer, intent(in)                     :: kold
      real(wp), dimension(neqn, 16), intent(in) :: phi
      integer, intent(in)                     :: ivc
      integer, dimension(10), intent(in)       :: iv
      integer, intent(in)                     :: kgi
      real(wp), dimension(11), intent(in)      :: gi
      real(wp), dimension(12), intent(in)      :: alpha
      real(wp), dimension(13), intent(in)      :: og
      real(wp), dimension(12), intent(in)      :: ow
      real(wp), dimension(neqn), intent(in)    :: oy
      real(wp), intent(in)                    :: ox

      !local variables:
      integer  :: i, iq, iw, j, jq, kp1, kp2, l, m
      real(wp), dimension(13) :: c, g, w
      real(wp) :: alp, gdi, gdif, gamma, h, hi, &
                  hmu, rmu, sigma, temp1, temp2, temp3, &
                  xi, xim1, xiq

      kp1 = kold + 1
      kp2 = kold + 2
      hi = xout - ox
      h = x - ox
      xi = hi/h
      xim1 = xi - 1.0_wp

      ! initialize w(*) for computing g(*)
      xiq = xi
      do iq = 1, kp1
         xiq = xi*xiq
         temp1 = iq*(iq + 1)
         w(iq) = xiq/temp1
      end do

      ! compute the double integral term gdi
      if (kold <= kgi) then
         gdi = gi(kold)
      else
         if (ivc > 0) then
            iw = iv(ivc)
            gdi = ow(iw)
            m = kold - iw + 3
         else
            gdi = 1.0_wp/temp1
            m = 2
         end if
         if (m <= kold) then
            do i = m, kold
               gdi = ow(kp2 - i) - alpha(i)*gdi
            end do
         end if
      end if

      ! compute g(*) and c(*)
      g(1) = xi
      g(2) = 0.5_wp*xi*xi
      c(1) = 1.0_wp
      c(2) = xi
      if (kold >= 2) then
         do i = 2, kold
            alp = alpha(i)
            gamma = 1.0_wp + xim1*alp
            l = kp2 - i
            do jq = 1, l
               w(jq) = gamma*w(jq) - alp*w(jq + 1)
            end do
            g(i + 1) = w(1)
            c(i + 1) = gamma*c(i)
         end do
      end if

      ! define interpolation parameters
      sigma = (w(2) - xim1*w(1))/gdi
      rmu = xim1*c(kp1)/gdi
      hmu = rmu/h

      ! interpolate for the solution -- yout
      ! and for the derivative of the solution -- ypout
      yout = 0.0_wp
      ypout = 0.0_wp
      do j = 1, kold
         i = kp2 - j
         gdif = og(i) - og(i - 1)
         temp2 = (g(i) - g(i - 1)) - sigma*gdif
         temp3 = (c(i) - c(i - 1)) + rmu*gdif
         do l = 1, neqn
            yout(l) = yout(l) + temp2*phi(l, i)
            ypout(l) = ypout(l) + temp3*phi(l, i)
         end do
      end do
      do l = 1, neqn
         yout(l) = ((1.0_wp - sigma)*oy(l) + sigma*y(l)) + &
                   h*(yout(l) + (g(1) - sigma*og(1))*phi(l, 1))
         ypout(l) = hmu*(oy(l) - y(l)) + &
                    (ypout(l) + (c(1) + rmu*og(1))*phi(l, 1))
      end do

   end subroutine dintp