 |
Celeritas 0.6.0-dev.115+3b60a5fd
|
Loading...
Searching...
No Matches
|
Celeritas 0.6.0-dev.115+3b60a5fd
|
Calculate the second derivatives of a cubic spline. More...
#include <SplineDerivCalculator.hh>
Public Types | |
| enum class | BoundaryCondition { natural = 0 , not_a_knot , geant , size_ } |
| Cubic spline interpolation boundary conditions. More... | |
Type aliases | |
| using | SpanConstReal = detail::SpanGridAccessor::SpanConstReal |
| using | Values = detail::XsGridAccessor::Values |
| using | VecReal = std::vector< real_type > |
Public Member Functions | |
| SplineDerivCalculator (BoundaryCondition) | |
| Construct from boundary conditions. | |
| VecReal | operator() (SpanConstReal, SpanConstReal) const |
| Calculate the second derivatives from spans. | |
| VecReal | operator() (XsGridData const &, Values const &) const |
| Calculate the second derivatives from grid data. | |
| template<class GA > | |
| auto | operator() (GA &&grid) const -> VecReal |
| Calculate the second derivatives. | |
| template<class GA > | |
| void | calc_initial_coeffs (GA const &grid, Real3 &tridiag, real_type &rhs) const |
| Calculate the coefficients for the first row using the boundary conditions. | |
| template<class GA > | |
| void | calc_final_coeffs (GA const &grid, Real3 &tridiag, real_type &rhs) const |
| Calculate the coefficients for the last row using the boundary conditions. | |
| template<class GA > | |
| void | calc_boundaries (GA const &grid, VecReal &deriv) const |
| Calculate the first and last values of the second derivative. | |
| template<class GA > | |
| auto | calc_geant_derivatives (GA const &grid) const -> VecReal |
| Calculate the second derivatives. | |
Calculate the second derivatives of a cubic spline.
See section 3.3: Cubic Spline Interpolation in [press-nr-1992] for a review of interpolating cubic splines and an algorithm for calculating the second derivatives.
Determining the polynomial coefficients $ a_0, a_1, a_2 \( and \$f a_3 \) of a cubic spline \( S(x) \) (see:
\[ h_{i - 1} S''_{i - 1} + 2 (h_{i - 1} + h_i) S''i + h_i S''_{i + 1} = 6 r_i, \]
where \( r_i = frac{\Delta y_i}{h_i} - frac{\Delta y_{i - 1}}{h_{i - 1}} \) and \( h_i = \Delta x_i = x_{i + 1} - x_i \).Specifying the boundary conditions gives the remaining two equations. Natural boundary conditions set \( S''_0 = S''_{n - 1} = 0 \), which leads to the following initial and final equations:
\begin{align} 2 (h_0 + h_1) S''_1 + h_1 S''_2 &= 6 r_1 // h_{n - 3} S''_{n - 3} + 2 (h_{n - 3} + h_{n - 2}) S''_{n - 2} &= 6 r_{n - 2}. \end{align}
The points \( x_0, x_1, \dots , x_{n - 1} \) where the spline changes from one cubic to the next are called knots. "Not-a-knot" boundary conditions require the third derivative \( S'''_i \) to be continous across the first and final interior knots, \( x_1 \) and \( x_{n - 2} \) (the name refers to the polynomials on the interval \( (x_0, x_1) \) and \( (x_1, x_2) \) being the same cubic, so \( x_1 \) is "not a knot"). This constraint gives the final two equations:
\begin{align} \frac{(h_0 + 2 h_1)(h_0 + h_1)}{h_1} S''_1 + \frac{h_1^2 - h_0^2}{h_1} S''_2 &= 6 r_1 // \frac{h_{n - 3}^2 - h_{n - 2}^2}{h_{n - 3}} S''_{n - 3} + \frac{(h_{n - 3} + h_{n - 2})(2 h_{n - 3} + h_{n - 2})}{h_{n - 3}} S''_{n - 2} &= 6 r_{n - 2} \end{align}
Once the system of equations has been solved for the second derivatives, the derivatives \( S''_0 \) and \( S''_{n - 1} \) can be recovered:
\begin{align} S''_0 &= \frac{(h_0 + h_1) S''_1 - h_0 S''_2}{h_1} \\ S''_{n - 1} &= \frac{(h_{n - 3} + h_{n - 2}) S''_{n - 2} - h_{n - 2} S''_{n - 3}}{h_{n - 3}} \end{align}
Calculate the second derivatives.
This is a hack to produce the same interpolation results as Geant4. The calculation here is identical to Geant4's G4PhysicsVector::ComputeSecDerivative1, which is based off the algorithm for calculating the second derivatives of a cubic spline in Numerical Recipes, modified for not-a-knot boundary conditions.
Note that here the coefficients are divided by \( h_i + h_{i + 1} \).
1.9.8