subroutine interpolate_r1(comm,f,t_want,y_want,yderiv_want)
!
! Part of rksuite_90 v1.0 (Aug 1994)
! software for initial value problems in ODEs
!
! Authors: R.W. Brankin (NAG Ltd., Oxford, England)
! I. Gladwell (Math Dept., SMU, Dallas, TX, USA)
! see main doc for contact details
!
real(kind=wp), intent(in) :: t_want !indep!
type(rk_comm_real_1d), intent(inout), target :: comm
real(kind=wp), dimension(:), intent(out), optional :: y_want !dep!
real(kind=wp), dimension(:), intent(out), optional :: yderiv_want !dep!
!
interface
function f(t,y)
use rksuite_90_prec, only:wp
real(kind=wp), intent(in) :: t !indep!
real(kind=wp), dimension(:), intent(in) :: y !dep!
real(kind=wp), dimension(size(y,1)) :: f !dep!
end function f
end interface
!
character(len=*),parameter :: srname="INTERPOLATE"
integer :: ier, jer, nrec, state, npcls
logical :: intrp_initialised
!
integer, parameter :: not_ready=-1, usable=-2, fatal=911, catastrophe=912, &
just_fine=1
logical, parameter :: ask=.true.
!
ier = just_fine; nrec = 0
!
body: do
!
! Is it permissible to call INTERPOLATE?
!
state = get_saved_state_r1("STEP_INTEGRATE",comm%save_states)
if (state==fatal) then
ier = catastrophe; nrec = 1; write (comm%rec,"(a)") &
" ** A catastrophic error has already been detected elsewhere."
exit body
end if
if (comm%use_range) then
if (get_saved_state_r1("RANGE_INTEGRATE",comm%save_states)/=usable) then
ier = fatal; nrec = 2; write (comm%rec,"(a/a)") &
" ** You have called INTERPOLATE after you specified to SETUP that you were",&
" ** going to use RANGE_INTEGRATE. This is not permitted."
exit body
end if
end if
if (state==not_ready) then
ier = fatal; nrec = 1; write (comm%rec,"(a)") &
" ** You have not called STEP_INTEGRATE, so you cannot use INTERPOLATE."
exit body
end if
if (state > just_fine) then
ier = fatal; nrec = 2; write (comm%rec,"(a/a)") &
" ** STEP_INTEGRATE has returned with a flag value greater than 1. You", &
" ** cannot call INTERPOLATE in this circumstance."
exit body
end if
!
! Check sizes of arrays
!
if (present(y_want)) then
if (any(shape(y_want)/=shape(comm%y))) then
ier = fatal; nrec = 3; write (comm%rec,"(a,i6,a/a,i6,a/a)") &
" ** The shape of the array Y_WANT is not consistent with the shape of the ", &
" ** dependent variables."
exit body
end if
end if
if (present(yderiv_want)) then
if (any(shape(yderiv_want)/=shape(comm%y))) then
ier = fatal; nrec = 3; write (comm%rec,"(a,i6,a/a,i6,a/a)") &
" ** The shape of the array YDERIV_WANT is not consistent with the shape of", &
" ** the dependent variables."
exit body
end if
end if
!
! Check METHOD is ok to interpolate with
!
if (comm%rk_method==3) then
ier = fatal; nrec = 5; write (comm%rec,"(a/a/a/a/a)") &
" ** You have been using STEP_INTEGRATE with METHOD = 'H' to integrate your",&
" ** equations. You have just called INTERPOLATE, but interpolation is not",&
" ** available for this METHOD. Either use METHOD = 'M', for which",&
" ** interpolation is available, or use RESET_T_END to make STEP_INTEGRATE",&
" ** step exactly to the points where you want output."
exit body
end if
!
! Get some workspace -
! can overwrite STAGES in METHOD 'L' since they're not requird for the
! interpolant
!
select case(comm%rk_method)
case(1)
if (.not. associated(comm%p)) comm%p => comm%stages(:,1:2)
npcls = 2
if (.not. associated(comm%ytemp)) comm%p => comm%stages(:,1:3)
case(2)
jer = 0
if (.not.associated(comm%xstage)) then
allocate(comm%xstage(size(comm%y,1)),stat=jer) !alloc!
end if
if (.not.associated(comm%ytemp)) then
allocate(comm%ytemp(size(comm%y,1)),stat=jer) !alloc!
end if
if (.not.associated(comm%p)) then
allocate(comm%p(size(comm%y,1),5),stat=jer) !alloc!
end if
npcls = 5
if (jer /= 0) then
ier = fatal; nrec = 1 ; write (comm%rec,"(a)") &
" ** Not enough storage available to create workspace required internally."
exit body
end if
end select
!
! Check data to see if interpolant has already been calculated for this
! step
!
intrp_initialised = get_saved_state_r1(srname,comm%save_states) /= usable
!
! Some initialization must be done before interpolation is possible.
!
if (.not.intrp_initialised) call form_intrp(f,comm%p)
!
! The actual evaluation of the interpolating polynomial and/or its first
! derivative is done in EVALUATE_INTRP.
!
call evaluate_intrp(comm%p,y_want,yderiv_want)
exit body
!
end do body
!
call rkmsg_r1(ier,srname,nrec,comm)
!
contains
!
subroutine form_intrp(f,p)
!
real(kind=wp), intent(out), dimension(:,:) :: p !dep!
!
interface
function f(t,y)
use rksuite_90_prec, only:wp
real(kind=wp), intent(in) :: t !indep!
real(kind=wp), dimension(:), intent(in) :: y !dep!
real(kind=wp), dimension(size(y,1)) :: f !dep!
end function f
end interface
!
real(kind=wp), dimension(:,:), pointer :: r !real!
real(kind=wp), dimension(:,:), pointer :: stages !dep!
real(kind=wp), dimension(:), pointer :: y, yp, y_old, yp_old !dep!
real(kind=wp), dimension(:), pointer :: xstage !dep!
!
stages => comm%stages
r => comm%r
y => comm%y
yp => comm%yp
y_old => comm%y_old
yp_old => comm%yp_old
xstage => comm%xstage
!
select case(comm%rk_method)
case(1)
!
! METHOD = 'L'. Use the cubic Hermite interpolant that is fully
! specified by the values and slopes at the two ends of the step.
!
p(:,2) = y - y_old
p(:,1) = comm%h_old*yp - p(:,2)
p(:,2) = p(:,1) - (p(:,2)-comm%h_old*yp_old)
p(:,1) = p(:,1) + p(:,2)
!
case(2)
!
! METHOD = 'M'.
!
if (.not.intrp_initialised) call extra_stages(f,comm%ytemp,comm%xstage)
!
! Form the coefficients of the interpolating polynomial in its shifted
! and scaled form. The transformation from the form in which the
! polynomial is derived can be somewhat ill-conditioned. The terms
! are grouped so as to minimize the errors of the transformation.
!
! Coefficient of SIGMA**6
p(:,5) = r(5,6)*stages(:,4) + &
((r(10,6)*xstage+r(8,6)*yp)+ &
(r(7,6)*stages(:,6)+r(6,6)*stages(:,5))) + &
((r(4,6)*stages(:,3)+r(9,6)*stages(:,7))+ &
(r(3,6)*stages(:,2)+r(11,6)*stages(:,1))+ &
r(1,6)*yp_old)
!
! Coefficient of SIGMA**5
p(:,4) = (r(10,5)*xstage+r(9,5)*stages(:,7)) + &
((r(7,5)*stages(:,6)+r(6,5)*stages(:,5))+ &
r(5,5)*stages(:,4)) + ((r(4,5)*stages(:,3)+ &
r(8,5)*yp)+(r(3,5)*stages(:,2)+r(11,5)* &
stages(:,1))+r(1,5)*yp_old)
!
! Coefficient of SIGMA**4
p(:,3) = ((r(4,4)*stages(:,3)+r(8,4)*yp)+ &
(r(7,4)*stages(:,6)+r(6,4)*stages(:,5))+ &
r(5,4)*stages(:,4)) + ((r(10,4)*xstage+ &
r(9,4)*stages(:,7))+(r(3,4)*stages(:,2)+ &
r(11,4)*stages(:,1))+r(1,4)*yp_old)
!
! Coefficient of SIGMA**3
p(:,2) = r(5,3)*stages(:,4) + r(6,3)*stages(:,5) + &
((r(3,3)*stages(:,2)+r(9,3)*stages(:,7))+ &
(r(10,3)*xstage+r(8,3)*yp)+r(1,3)* &
yp_old)+((r(4,3)*stages(:,3)+r(11,3)* &
stages(:,1))+r(7,3)*stages(:,6))
!
! Coefficient of SIGMA**2
p(:,1) = r(5,2)*stages(:,4) + ((r(6,2)*stages(:,5)+ &
r(8,2)*yp)+r(1,2)*yp_old) + &
((r(3,2)*stages(:,2)+r(9,2)*stages(:,7))+ &
r(10,2)*xstage) + ((r(4,2)*stages(:,3)+ &
r(11,2)*stages(:,1))+r(7,2)*stages(:,6))
!
! Scale all the coefficients by the step size.
p(:,:) = comm%h_old*p(:,:)
!
end select
!
end subroutine form_intrp
subroutine evaluate_intrp(p,y_want,yderiv_want)
!
real(kind=wp), dimension(:), optional, intent(out) :: y_want !dep!
real(kind=wp), dimension(:), optional, intent(out) :: yderiv_want !dep!
real(kind=wp), dimension(:,:), intent(in) :: p !dep!
!
real :: sigma
integer :: i
!
sigma = (t_want-comm%t)/comm%h_old
!
if (present(y_want)) then
y_want = p(:,comm%intrp_degree-1)*sigma
do i = comm%intrp_degree - 2, 1, -1
y_want = (y_want+p(:,i))*sigma
end do
y_want = (y_want+comm%h_old*comm%yp)*sigma + comm%y
end if
!
if (present(yderiv_want)) then
yderiv_want = comm%intrp_degree*p(:,comm%intrp_degree-1)*sigma
do i = comm%intrp_degree - 1, 2, -1
yderiv_want = (yderiv_want+i*p(:,i-1))*sigma
end do
yderiv_want = (yderiv_want+comm%h_old*comm%yp)/comm%h_old
end if
!
end subroutine evaluate_intrp
subroutine extra_stages(f,ytemp,xstage)
!
real(kind=wp), dimension(:), intent(out) :: ytemp, xstage !dep!
!
interface
function f(t,y)
use rksuite_90_prec, only:wp
real(kind=wp), intent(in) :: t !indep!
real(kind=wp), dimension(:), intent(in) :: y !dep!
real(kind=wp), dimension(size(y,1)) :: f !dep!
end function f
end interface
!
real(kind=wp), dimension(:,:), pointer :: stages !dep!
real(kind=wp), dimension(:), pointer :: yp, y_old, yp_old !dep!
!
real(kind=wp), dimension(:,:), pointer :: a !real!
real(kind=wp), dimension(:), pointer :: c !real!
real(kind=wp), pointer :: h_old, t_old !indep!
!
integer :: i, j
!
a => comm%a
stages => comm%stages
c => comm%c
yp => comm%yp
y_old => comm%y_old
yp_old => comm%yp_old
h_old => comm%h_old
t_old => comm%t_old
!
! Compute the extra stages needed for interpolation using the facts that
! 1. Stage 1 is YP_OLD.
! 2. Stage i (i>1) is stored in STAGES(...,i-1).
! 3. This pair is FSAL, i.e. STAGES(...,7)=YP, which frees
! up STAGES(...,7) for use by stage 9.
! 4. XSTAGE is used for stage 10.
! 5. The coefficient of stage 2 in the interpolant is always 0, so
! STAGES(...,1) is used for stage 11.
!
do i = 9, 11
do j = 1, i-1
select case (j)
case(1); ytemp = a(i,1)*yp_old
! could have used matmul here but that prevents increasing rank of dep-var
case(2:7);ytemp = ytemp + a(i,j)*stages(:,j-1)
case(8); ytemp = ytemp + a(i,j)*yp
case(9); ytemp = ytemp + a(i,j)*stages(:,7)
case(10); ytemp = ytemp + a(i,j)*xstage
end select
end do
ytemp = y_old + h_old*ytemp
select case(i)
case(9)
stages(:,7) = f(t_old+c(i)*h_old,ytemp)
case(10)
xstage = f(t_old+c(i)*h_old,ytemp)
case(11)
stages(:,1) = f(t_old+c(i)*h_old,ytemp)
end select
comm%f_count = comm%f_count + 1
end do
!
end subroutine extra_stages
end subroutine interpolate_r1