HQSVAL

This rank-1 array function evaluates a Hermite quintic spline or one of its derivatives for an array of input points. In particular, it computes solution values for the Feynman-Kac PDE handled by routine FEYNMAN_KAC.

Function Return Value

HQSVAL — Rank-1 array of extent  size(XVAL) containing the values or derivatives of the Hermite quintic spline at the points given in array XVAL. (Output)

Required Arguments

XVAL — Rank-1 array containing the points at which the Hermite quintic spline is to be evaluated. (Input)
Let NXVAL = size(XVAL). The points in XVAL must lie within the range of array BREAK, i.e. BREAK(1) ≤ XVAL(I) ≤ BREAK(NXVAL), I=1, …, NXVAL.

BREAK — Rank-1 array containing the breakpoints for the Hermite quintic spline representation. (Input)
When applied to the output from routine FEYNMAN_KAC, array BREAK is identical to array XGRID.
Let NBREAK = size(BREAK). NBREAK must be at least 2 and the values in BREAK must be in strictly increasing order.

COEFFS — Rank-1 array of size 3 * NBREAK containing the coefficients of the Hermite quintic spline representation. (Input)
When applied to the output arrays Y or YPRIME from routine FEYNMAN_KAC, array COEFFS is identical to one of the columns of arrays Y or YPRIME, respectively.

Optional Argument

IDERIV — Order of the derivative to be evaluated. (Input)
It is required that IDERIV = 0, 1, 2 or 3. Use 0 for the function itself, 1 for the first derivative, etc.
Default: IDERIV = 0.

FORTRAN 90 Interface

Generic: HQSVAL (XVAL, BREAK, COEFFS [, …])

Specific: The specific interface names are S_HQSVAL and D_HQSVAL.

Description

The Hermite quintic spline interpolation is done over the composite interval [xmin,xmax], where BREAK(I) = xi are given by (xmin =) x1 < x2 < … < xm (= xmax).

The Hermite quintic spline function is constructed using three primary functions, defined by

 

For each

 

the spline is locally defined by

 

where

 

are the values of a given twice continuously differentiable function f and its first two derivatives at the breakpoints. The approximating function H(x) is twice continuously differentiable on [xmin,xmax], whereas the third derivative is in general only continuous within the interior of the intervals [xi,xi+1]. From the local representation of H(x) it follows that

 

The spline coefficients yi, i = 1, .…, 3m are stored as successive triplets in array COEFFS. For a given w ∈ [xmin,xmax], function HQSVAL uses the information in COEFFS together with the values of b0, b1, b2 and its derivatives at w to compute H(d)(w), d = 0, … 3 using the local representation on the particular subinterval containing w.

When using arrays XGRID and Y(:,I) from routine FEYNMAN_KAC as input arrays BREAK and COEFFS, function HQSVAL allows for computation of approximate solutions to the Feynman-Kac PDE for IDERIV=0, 1, 2, 3, respectively. The solution values are evaluated at the array of points (XVAL(:),TGRID(I)) for I=1, …,size(TGRID) and (XVAL(:),0) for I=0 . Similarly, using arrays XGRID and YPRIME(:,I) allows for computation of approximate solutions to the Feynman-Kac PDE.

Example: Exact Interpolation with Hermite Quintic

Consider function f(x) = x5, a polynomial of degree 5, on the interval [‑1,1] with breakpoints ±1. Then, the end derivative values are

 

and

 

Since the Hermite quintic interpolates all polynomials up to degree 5 exactly, the spline interpolation on [‑1,1] must agree with the exact function values up to rounding errors.

 

use hqsval_int

use umach_int

 

implicit none

 

integer :: i, nout

real(kind(1e0)) :: break(2), xval(7), coeffs(6), interpolant(7)

 

! Define arrays

break = (/ -1.0, 1.0 /)

xval = (/ -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75 /)

coeffs = (/ -1.0, 5.0, -20.0, 1.0, 5.0, 20.0 /)

 

! Compute interpolant

interpolant = hqsval(xval, break, coeffs)

 

call umach(2, nout)

! Compare interpolant with exact function.

write(nout,'(A6,A10,A15,A10)')'x', 'F(x)', 'Interpolant', 'Error'

write(nout,'(f8.3,f9.3,f11.3,f15.7)') (xval(i), F(xval(i)), &

interpolant(i), abs(F(xval(i))-interpolant(i)), &

i=1,7)

 

contains

function F(x)

implicit none

real(kind(1e0)) :: F, x

 

F = x**5

return

end function F

end

 

Output

 

x F(x) Interpolant Error

-0.750 -0.237 -0.237 0.0000000

-0.500 -0.031 -0.031 0.0000000

-0.250 -0.001 -0.001 0.0000000

0.000 0.000 0.000 0.0000000

0.250 0.001 0.001 0.0000000

0.500 0.031 0.031 0.0000000

0.750 0.237 0.237 0.0000000