feynmanKacEvaluate

Computes the value of a Hermite quintic spline or the value of one of its derivatives. In particular, computes solutions to the Feynman-Kac PDE handled by function feynmanKac.

Synopsis

feynmanKacEvaluate (breakpoints, w, coef)

Required Arguments

float breakpoints[] (Input)

Array of length m containing the breakpoints for the Hermite quintic spline interpolation. The breakpoints must be in strictly increasing order. When applied to feynmanKac, breakpoints[] is identical to array xgrid[].

float w[] (Input)

Vector of length nw containing the evaluation points for the spline. It is required that breakpoints[0] $\leq$ w[i] $\leq$ breakpoints[m-1] for i=0,…, nw-1.

float coef[] (Input)

Vector of length 3*m containing the coefficients of the Hermite quintic spline.

When applied to feynmanKac, this vector is one of the rows of output arrays y or yPrime related to the spline coefficients at time points t=tgrid[j], j=1, …, ntgrid.

Return Value

An array of length nw containing the values of the Hermite quintic spline or one of its derivatives at the evaluation points in array w[]. If no values can be computed, then None is returned.

Optional Arguments

deriv, int (Input)

Let d = deriv and let $H(w)$ be the given Hermite quintic spline. Then, this option produces the d-th derivative of $H(w)$ at w, $H^d(w)$. It is required that deriv = 0,1,2 or 3.

Default: deriv = 0.

Description

The Hermite quintic spline interpolation is done over the composite interval $\left[x_{\text{min}}, x_{\text{max}}\right]$, where

breakpoints[i-1] = $x_i$ are given by $\left( x_{\min} = \right) x_1 < x_2 < \ldots < x_m \left( = x_{\max} \right)$ .

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

$$\begin{split}\begin{aligned} b_0(z) &= -6z^5 + 15z^4 - 10z^3 + 1 = (1-z)^3\left(6z^2 + 3z + 1\right), \\ b_1(z) &= -3z^5 + 8z^4 - 6z^3 + z = (1-z)^3z(3z+1), \\ b_2(z) &= \tfrac{1}{2}\left(-z^5+3z^4-3z^3+z^2\right) = \tfrac{1}{2}(1-z)^3z^2. \end{aligned}\end{split}$$

For each

$$x \in \left[x_i,x_{i+1}\right], h_i = x_{i+1} - x_i, z_i = \left(x-x_i\right) / h_i, i = 1, \ldots, m-1,$$

the spline is locally defined by

$$\begin{split}\begin{array}{c} H(x) = y_{3i-2}b_0(z) + y_{3i+1}b_0(1-z) + h_iy_{3i-1}b_1(z) \\ -h_iy_{3i+2}b_1(1-z) + h_i^2y_{3i}b_2(z) + h_i^2y_{3i+3}b_2(1-z), \\ \end{array}\end{split}$$

where

$$\begin{split}\begin{aligned} y_{3i-2} = f\left(x_i\right), y_{3i-1} = \left({\partial}f / {\partial}x\right)\left(x_i\right) &= f'\left(x_i\right), y_{3i} = \left(\partial^2f / {\partial}x^2\right)\left(x_i\right) \\ &= f''\left(x_i\right), i = 1, \ldots, m - 1 \end{aligned}\end{split}$$

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 $\left[x_{\text{min}}, x_{\text{max}}\right]$ , whereas the third derivative is in general only continuous within the interior of the intervals $\left[ x_i, x_{i+1} \right]$ . From the local representation of $H(x)$ it follows that

$$H\left(x_i\right) = f\left(x_i\right) = y_{3i-2}, H'\left(x_i\right) = f'\left(x_i\right) = y_{3i-1}, H''\left(x_i\right) = y_{3i}, i = 1, \ldots, m$$

The spline coefficients $y_i, i = 1, \ldots, 3m$ are stored as successive triplets in array coef[]. For a given $w \in \left[ x_{\min}, x_{\max} \right]$, function feynmanKacEvaluate uses the information in coef[] together with the values of $b_0, b_1, b_2$ and its derivatives at $w$ to compute $H^{(d)} (w), d = 0, \ldots, 3$ using the local representation on the particular subinterval containing $w$.

Example

Consider function $f(x) = x^5$, a polynomial of degree 5, on the interval $[-1, 1]$ with breakpoints $\pm 1$. Then, the end derivative values are

$$y_1 = f(-1) = -1, y_2 = f'(-1) = 5, y_3 = f''(-1) = -20$$

and

$$y_4 = f(1) = 1, y_5 = f'(1) = 5, y_6 = f''(1) = 20$$

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

from __future__ import print_function
from numpy import *
from pyimsl.math.feynmanKacEvaluate import feynmanKacEvaluate

# Define function


def func(x):
    return x ** 5.0


nw = 7
m = 2
breakpoints = [-1.0, 1.0]
w = [-0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75]
coef = [-1.0, 5.0, -20.0, 1.0, 5.0, 20.0]

result = feynmanKacEvaluate(breakpoints, w, coef)

print("     x          F(x)   Interpolant   Error\n")
for i in range(0, 7):
    print("  %6.3f  %10.3f  %10.3f  %10.7f" %
          (w[i], func(w[i]), result[i], fabs(func(w[i]) - result[i])))

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