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