Piecewise Cubic Hermite Interpolation to Monotone data.
Sets derivatives needed to determine a monotone piecewise cubic Hermite interpolant to the data given in X and F.
Default boundary conditions are provided which are compatible with monotonicity. (See dpchic if user control of boundary conditions is desired.)
If the data are only piecewise monotonic, the interpolant will have an extremum at each point where monotonicity switches direction. (See dpchic if user control is desired in such cases.)
To facilitate two-dimensional applications, includes an increment between successive values of the F- and D-arrays.
The resulting piecewise cubic Hermite function may be evaluated by dpchfe or dpchfd.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n |
number of data points. (Error return if N<2). If N=2, simply does linear interpolation. |
||
| real(kind=wp), | intent(in) | :: | x(*) |
array of independent variable values. The
elements of X must be strictly increasing:
|
||
| real(kind=wp), | intent(in) | :: | f(incfd,*) |
array of dependent variable values to be interpolated. F(1+(I-1)*INCFD) is value corresponding to X(I). DPCHIM is designed for monotonic data, but it will work for any F-array. It will force extrema at points where monotonicity switches direction. If some other treatment of switch points is desired, dpchic should be used instead. |
||
| real(kind=wp), | intent(out) | :: | d(incfd,*) |
array of derivative values at the data
points. If the data are monotonic, these values will
determine a monotone cubic Hermite function.
The value corresponding to X(I) is stored in
|
||
| integer, | intent(in) | :: | incfd |
increment between successive values in F and D. This argument is provided primarily for 2-D applications. (Error return if INCFD<1). |
||
| integer, | intent(out) | :: | ierr |
error flag.
NOTE: The above errors are checked in the order listed, and following arguments have NOT been validated. |
subroutine dpchim (n, x, f, d, incfd, ierr) integer,intent(in) :: n !! number of data points. (Error return if N<2). !! If N=2, simply does linear interpolation. integer,intent(in) :: incfd !! increment between successive values in F and D. !! This argument is provided primarily for 2-D applications. !! (Error return if INCFD<1). integer,intent(out) :: ierr !! error flag. !! !! * Normal return: !! * IERR = 0 (no errors). !! !! * Warning error: !! * IERR>0 means that IERR switches in the direction !! of monotonicity were detected. !! !! * "Recoverable" errors (The D-array has not been changed in any of these cases.): !! * IERR = -1 if N<2 . !! * IERR = -2 if INCFD<1 . !! * IERR = -3 if the X-array is not strictly increasing. !! !! NOTE: The above errors are checked in the order listed, !! and following arguments have **NOT** been validated. real(wp),intent(in) :: x(*) !! array of independent variable values. The !! elements of X must be strictly increasing: !! `X(I-1) < X(I), I = 2(1)N`. !! (Error return if not.) real(wp),intent(in) :: f(incfd,*) !! array of dependent variable values to be !! interpolated. F(1+(I-1)*INCFD) is value corresponding to !! X(I). DPCHIM is designed for monotonic data, but it will !! work for any F-array. It will force extrema at points where !! monotonicity switches direction. If some other treatment of !! switch points is desired, [[dpchic]] should be used instead. real(wp),intent(out) :: d(incfd,*) !! array of derivative values at the data !! points. If the data are monotonic, these values will !! determine a monotone cubic Hermite function. !! The value corresponding to X(I) is stored in !! `D(1+(I-1)*INCFD), I=1(1)N`. !! No other entries in D are changed. integer :: i, nless1 real(wp) :: del1, del2, dmax, dmin, drat1, drat2, dsave, & h1, h2, hsum, hsumt3, w1, w2, tmp ! validity-check arguments. if ( n<2 ) then ierr = -1 call xermsg ('PCHIP', 'dpchim', 'number of data points less than two', ierr, 1) return end if if ( incfd<1 ) then ierr = -2 call xermsg ('PCHIP', 'dpchim', 'increment less than one', ierr, 1) return end if do i = 2, n if ( x(i)<=x(i-1) ) then ! x-array not strictly increasing. ierr = -3 call xermsg ('PCHIP', 'dpchim', 'x-array not strictly increasing', ierr, 1) end if end do ! function definition is ok, go on. ierr = 0 nless1 = n - 1 h1 = x(2) - x(1) del1 = (f(1,2) - f(1,1))/h1 dsave = del1 ! special case n=2 -- use linear interpolation. if (nless1 <= 1) then d(1,1) = del1 d(1,n) = del1 return end if ! normal case (n >= 3). h2 = x(3) - x(2) del2 = (f(1,3) - f(1,2))/h2 ! set d(1) via non-centered three-point formula, adjusted to be ! shape-preserving. hsum = h1 + h2 w1 = (h1 + hsum)/hsum w2 = -h1/hsum d(1,1) = w1*del1 + w2*del2 if ( dpchst(d(1,1),del1) <= zero) then d(1,1) = zero else if ( dpchst(del1,del2) < zero) then ! need do this check only if monotonicity switches. dmax = three*del1 if (abs(d(1,1)) > abs(dmax)) d(1,1) = dmax endif ! loop through interior points. do i = 2, nless1 if (i /= 2) then h1 = h2 h2 = x(i+1) - x(i) hsum = h1 + h2 del1 = del2 del2 = (f(1,i+1) - f(1,i))/h2 end if ! set d(i)=0 unless data are strictly monotonic. d(1,i) = zero tmp = dpchst(del1,del2) if (tmp < zero) then ierr = ierr + 1 dsave = del2 elseif (tmp == zero) then ! count number of changes in direction of monotonicity. if (del2 /= zero) then if ( dpchst(dsave,del2) < zero) ierr = ierr + 1 dsave = del2 end if else ! use brodlie modification of butland formula. hsumt3 = hsum+hsum+hsum w1 = (hsum + h1)/hsumt3 w2 = (hsum + h2)/hsumt3 dmax = max( abs(del1), abs(del2) ) dmin = min( abs(del1), abs(del2) ) drat1 = del1/dmax drat2 = del2/dmax d(1,i) = dmin/(w1*drat1 + w2*drat2) end if end do ! set d(n) via non-centered three-point formula, adjusted to be ! shape-preserving. w1 = -h2/hsum w2 = (h2 + hsum)/hsum d(1,n) = w1*del1 + w2*del2 if ( dpchst(d(1,n),del2) <= zero) then d(1,n) = zero else if ( dpchst(del1,del2) < zero) then ! need do this check only if monotonicity switches. dmax = three*del2 if (abs(d(1,n)) > abs(dmax)) d(1,n) = dmax endif end subroutine dpchim