dvset is called by dvstep and sets coefficients for use there.
for each order nq, the coefficients in el are calculated by use of
the generating polynomial lambda(x), with coefficients el(i).
for the backward differentiation formulas,
for the adams formulas,
where c is a normalization constant.
in both cases, is defined by
in addition to variables described previously, communication with dvset uses the following class variables:
tau = a vector of length 13 containing the past nq values
of h.el = a vector of length 13 in which vset stores the
coefficients for the corrector formula.tq = a vector of length 5 in which vset stores constants
used for the convergence test, the error test, and the
selection of h at a new order.meth = the basic method indicator.nq = the current order.l = nq + 1, the length of the vector stored in el, and
the number of columns of the yh array being used.nqwait = a counter controlling the frequency of order changes.
an order change is about to be considered if nqwait = 1.| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(dvode_t), | intent(inout) | :: | me |
subroutine dvset(me) class(dvode_t),intent(inout) :: me real(wp) :: ahatn0 , alph0 , cnqm1 , csum , elp , & em(13) , em0 , floti , flotl , flotnq , hsum , & rxi , rxis , s , t1 , t2 , t3 , t4 , t5 , & t6 , xi integer :: i , iback , j , jp1 , nqm1 , nqm2 real(wp),parameter :: cortes = 0.1_wp real(wp),parameter :: one = 1.0_wp real(wp),parameter :: six = 6.0_wp real(wp),parameter :: two = 2.0_wp flotl = real(me%dat%l,wp) nqm1 = me%dat%nq - 1 nqm2 = me%dat%nq - 2 if ( me%dat%meth==2 ) then ! set coefficients for bdf methods. ------------------------------------ do i = 3 , me%dat%l me%dat%el(i) = zero enddo me%dat%el(1) = one me%dat%el(2) = one alph0 = -one ahatn0 = -one hsum = me%dat%h rxi = one rxis = one if ( me%dat%nq/=1 ) then do j = 1 , nqm2 ! in el, construct coefficients of (1+x/xi(1))*...*(1+x/xi(j+1)). ------ hsum = hsum + me%dat%tau(j) rxi = me%dat%h/hsum jp1 = j + 1 alph0 = alph0 - one/real(jp1,wp) do iback = 1 , jp1 i = (j+3) - iback me%dat%el(i) = me%dat%el(i) + me%dat%el(i-1)*rxi enddo enddo alph0 = alph0 - one/real(me%dat%nq,wp) rxis = -me%dat%el(2) - alph0 hsum = hsum + me%dat%tau(nqm1) rxi = me%dat%h/hsum ahatn0 = -me%dat%el(2) - rxi do iback = 1 , me%dat%nq i = (me%dat%nq+2) - iback me%dat%el(i) = me%dat%el(i) + me%dat%el(i-1)*rxis enddo endif t1 = one - ahatn0 + alph0 t2 = one + real(me%dat%nq,wp)*t1 me%dat%tq(2) = abs(alph0*t2/t1) me%dat%tq(5) = abs(t2/(me%dat%el(me%dat%l)*rxi/rxis)) if ( me%dat%nqwait==1 ) then cnqm1 = rxis/me%dat%el(me%dat%l) t3 = alph0 + one/real(me%dat%nq,wp) t4 = ahatn0 + rxi elp = t3/(one-t4+t3) me%dat%tq(1) = abs(elp/cnqm1) hsum = hsum + me%dat%tau(me%dat%nq) rxi = me%dat%h/hsum t5 = alph0 - one/real(me%dat%nq+1,wp) t6 = ahatn0 - rxi elp = t2/(one-t6+t5) me%dat%tq(3) = abs(elp*rxi*(flotl+one)*t5) endif ! set coefficients for adams methods. ---------------------------------- else if ( me%dat%nq/=1 ) then hsum = me%dat%h em(1) = one flotnq = flotl - one do i = 2 , me%dat%l em(i) = zero enddo do j = 1 , nqm1 if ( (j==nqm1) .and. (me%dat%nqwait==1) ) then s = one csum = zero do i = 1 , nqm1 csum = csum + s*em(i)/real(i+1,wp) s = -s enddo me%dat%tq(1) = em(nqm1)/(flotnq*csum) endif rxi = me%dat%h/hsum do iback = 1 , j i = (j+2) - iback em(i) = em(i) + em(i-1)*rxi enddo hsum = hsum + me%dat%tau(j) enddo ! compute integral from -1 to 0 of polynomial and of x times it. ------- s = one em0 = zero csum = zero do i = 1 , me%dat%nq floti = real(i,wp) em0 = em0 + s*em(i)/floti csum = csum + s*em(i)/(floti+one) s = -s enddo ! in el, form coefficients of normalized integrated polynomial. -------- s = one/em0 me%dat%el(1) = one do i = 1 , me%dat%nq me%dat%el(i+1) = s*em(i)/real(i,wp) enddo xi = hsum/me%dat%h me%dat%tq(2) = xi*em0/csum me%dat%tq(5) = xi/me%dat%el(me%dat%l) if ( me%dat%nqwait==1 ) then ! for higher order control constant, multiply polynomial by 1+x/xi(q). - rxi = one/xi do iback = 1 , me%dat%nq i = (me%dat%l+1) - iback em(i) = em(i) + em(i-1)*rxi enddo ! compute integral of polynomial. -------------------------------------- s = one csum = zero do i = 1 , me%dat%l csum = csum + s*em(i)/real(i+1,wp) s = -s enddo me%dat%tq(3) = flotl*em0/csum endif else me%dat%el(1) = one me%dat%el(2) = one me%dat%tq(1) = one me%dat%tq(2) = two me%dat%tq(3) = six*me%dat%tq(2) me%dat%tq(5) = one endif me%dat%tq(4) = cortes*me%dat%tq(2) end subroutine dvset