Cubic spline with a not-a-knot boundary conditions. More...
#include <spline.hpp>
Inherits sirius::Radial_grid< double >.
Public Member Functions | |
| Spline () | |
| Default constructor. More... | |
| Spline (Radial_grid< U > const &radial_grid__) | |
| Constructor of a new empty spline. More... | |
| Spline (Radial_grid< U > const &radial_grid__, std::function< T(U)> f__) | |
| Constructor of a spline from a function. More... | |
| Spline (Radial_grid< U > const &radial_grid__, std::vector< T > const &y__) | |
| Constructor of a spline from a list of values. More... | |
| Spline (Spline< T, U > &&src__)=default | |
| Spline< T, U > & | operator= (Spline< T, U > &&src__)=default |
| Spline< T, U > & | operator= (std::function< T(U)> f__) |
| Spline< T, U > & | operator= (std::vector< T > const &y__) |
| T & | operator() (const int i) |
| Get the reference to a value at the point x[i]. More... | |
| T const & | operator() (const int i) const |
| Get value at the point x[i]. More... | |
| T | operator() (const int i, U dx) const |
| Get value at the point x[i] + dx. More... | |
| T | at_point (U x) const |
| Compute value at any point. More... | |
| T | deriv (const int dm, const int i, const U dx) const |
| T | deriv (int dm, int i) const |
| Spline< T, U > & | interpolate () |
| T | integrate (std::vector< T > &g__, int m__) const |
| Integrate spline with r^m prefactor. More... | |
| T | integrate (int m__) const |
| Integrate spline with r^m prefactor. More... | |
| void | scale (double a__) |
| std::array< T, 4 > | coeffs (int i__) const |
| sddk::mdarray< T, 2 > const & | coeffs () const |
| std::vector< T > | values () const |
 Public Member Functions inherited from sirius::Radial_grid< double > | |
| Radial_grid () | |
| Constructor of the empty object. More... | |
| Radial_grid (int num_points__) | |
| Constructor. More... | |
| Radial_grid (Radial_grid< double > &&src__)=default | |
| Radial_grid< double > & | operator= (Radial_grid< double > &&src__)=default |
| int | index_of (double x__) const |
| int | num_points () const |
| Number of grid points. More... | |
| double | first () const |
| First point of the grid. More... | |
| double | last () const |
| Last point of the grid. More... | |
| double | operator[] (const int i) const |
| Return \( x_{i} \). More... | |
| double | x (const int i) const |
| Return \( x_{i} \). More... | |
| sddk::mdarray< double, 1 > const & | x () const |
| double | dx (const int i) const |
| Return \( dx_{i} \). More... | |
| sddk::mdarray< double, 1 > const & | dx () const |
| double | x_inv (const int i) const |
| Return \( x_{i}^{-1} \). More... | |
| std::string const & | name () const |
| void | copy_to_device () |
| Radial_grid< double > | segment (int num_points__) const |
| std::vector< real_type< double > > | values () const |
| uint64_t | hash () const |
Private Member Functions | |
| Spline (Spline< T, U > const &src__)=delete | |
| Spline< T, U > & | operator= (Spline< T, U > const &src__)=delete |
| int | solve (T *dl, T *d, T *du, T *b, int n) |
| Solver tridiagonal system of linear equaitons. More... | |
| void | init_grid (Radial_grid< U > const &radial_grid__) |
| Init the underlying radial grid. More... | |
Private Attributes | |
| sddk::mdarray< T, 2 > | coeffs_ |
| Array of spline coefficients. More... | |
Additional Inherited Members | |
| Protected Member Functions inherited from sirius::Radial_grid< double > | |
| void | init () |
| Initialize the grid. More... | |
| Radial_grid (Radial_grid< double > const &src__)=delete | |
| Radial_grid & | operator= (Radial_grid< double > const &src__)=delete |
| Protected Attributes inherited from sirius::Radial_grid< double > | |
| sddk::mdarray< double, 1 > | x_ |
| Radial grid points. More... | |
| sddk::mdarray< double, 1 > | x_inv_ |
| Inverse values of radial grid points. More... | |
| sddk::mdarray< double, 1 > | dx_ |
| Radial grid points difference. More... | |
| std::string | name_ |
| Name of the grid type. More... | |
Detailed Description
class sirius::Spline< T, U >
Cubic spline with a not-a-knot boundary conditions.
The following convention for spline coefficients is used: for \( x \) in \( [x_i, x_{i+1}] \) the value of the spline is equal to \( a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3 \).
Suppose we have \( n \) value points \( y_i = f(x_i) \) and \( n - 1 \) segments:
\[ S_i(x) = y_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3 \]
Segment derivatives:
\[ \begin{eqnarray} S_i'(x) &=& b_i + 2c_i(x - x_i) + 3d_i(x - x_i)^2 \\ S_i''(x) &=& 2c_i + 6d_i(x-x_i) \\ S_i'''(x) &=& 6d_i \end{eqnarray} \]
The following substitutions are made:
\[ m_i = 2 c_i \\ h_i = x_{i+1} - x_i \\ y_i' = \frac{y_{i+1} - y_i}{h_i} \]
Now we can equate the derivatives at the end points of segments. From the 3rd derivative we get
\[ d_i = \frac{1}{6h_i}(m_{i+1} - m_i) \]
From the 1st derivative we get
\[ b_i = y_i' - \frac{h_i}{2} m_i - \frac{h_i}{6}(m_{i+1} - m_i) \]
Using 2nd derivative condition we get
\[ h_i m_i + 2(h_{i} + h_{i+1})m_{i+1} + h_{i+1}m_{i+2} = 6(y_{i+1}' - y_{i}') \]
So far we got \( n - 3 \) equations for \( n - 1 \) coefficients \( m_i \). We need two extra conditions. Not-a-knot boundary condition (counting segments and points from 1):
\[ S_1'''(x_2) = S_2'''(x_2) \longrightarrow d_1 = d_2 \\ S_{n-2}'''(x_{n-1}) = S_{n-1}'''(x_{n-1}) \longrightarrow d_{n-2} = d_{n-1} \]
Definition at line 75 of file spline.hpp.
Constructor & Destructor Documentation
◆ Spline() [1/4]
|
inline |
Default constructor.
Definition at line 137 of file spline.hpp.
◆ Spline() [2/4]
|
inline |
Constructor of a new empty spline.
Definition at line 141 of file spline.hpp.
◆ Spline() [3/4]
|
inline |
Constructor of a spline from a function.
Definition at line 148 of file spline.hpp.
◆ Spline() [4/4]
|
inline |
Constructor of a spline from a list of values.
Definition at line 158 of file spline.hpp.
Member Function Documentation
◆ solve()
|
inlineprivate |
Solver tridiagonal system of linear equaitons.
Definition at line 85 of file spline.hpp.
◆ init_grid()
|
inlineprivate |
Init the underlying radial grid.
Definition at line 127 of file spline.hpp.
◆ operator=() [1/2]
|
inline |
Definition at line 175 of file spline.hpp.
◆ operator=() [2/2]
|
inline |
Definition at line 183 of file spline.hpp.
◆ operator()() [1/3]
|
inline |
Get the reference to a value at the point x[i].
Definition at line 195 of file spline.hpp.
◆ operator()() [2/3]
|
inline |
Get value at the point x[i].
Definition at line 201 of file spline.hpp.
◆ operator()() [3/3]
|
inline |
Get value at the point x[i] + dx.
Definition at line 207 of file spline.hpp.
◆ at_point()
|
inline |
Compute value at any point.
Definition at line 216 of file spline.hpp.
◆ deriv() [1/2]
|
inline |
Definition at line 231 of file spline.hpp.
◆ deriv() [2/2]
|
inline |
Definition at line 263 of file spline.hpp.
◆ interpolate()
|
inline |
- this should be boundary conditions for natural spline (with zero second derivatives at boundaries) */
Definition at line 275 of file spline.hpp.
◆ integrate() [1/2]
|
inline |
Integrate spline with r^m prefactor.
Derivation for r^2 prefactor is based on the following Mathematica notebook:
In[26]:= result =
FullSimplify[
Integrate[
x^(2)*(a0 + a1*(x - x0) + a2*(x - x0)^2 + a3*(x - x0)^3), {x, x0,
x1}],
Assumptions -> {Element[{x0, x1}, Reals], x1 > x0 > 0}]
Out[26]= 1/60 (20 a0 (-x0^3 + x1^3) +
5 a1 (x0^4 - 4 x0 x1^3 + 3 x1^4) + (x0 -
x1)^3 (-2 a2 (x0^2 + 3 x0 x1 + 6 x1^2) +
a3 (x0 - x1) (x0^2 + 4 x0 x1 + 10 x1^2)))
In[27]:= r = Expand[result] /. {x1 -> x0 + dx}
Out[27]= -((a0 x0^3)/3) + (a1 x0^4)/12 - (a2 x0^5)/30 + (
a3 x0^6)/60 + 1/3 a0 (dx + x0)^3 - 1/3 a1 x0 (dx + x0)^3 +
1/3 a2 x0^2 (dx + x0)^3 - 1/3 a3 x0^3 (dx + x0)^3 +
1/4 a1 (dx + x0)^4 - 1/2 a2 x0 (dx + x0)^4 +
3/4 a3 x0^2 (dx + x0)^4 + 1/5 a2 (dx + x0)^5 -
3/5 a3 x0 (dx + x0)^5 + 1/6 a3 (dx + x0)^6
In[34]:= Collect[r, dx, Simplify]
Out[34]= (a3 dx^6)/6 + a0 dx x0^2 + 1/2 dx^2 x0 (2 a0 + a1 x0) +
1/5 dx^5 (a2 + 2 a3 x0) + 1/3 dx^3 (a0 + x0 (2 a1 + a2 x0)) +
1/4 dx^4 (a1 + x0 (2 a2 + a3 x0))
In[28]:= r1 = Collect[r/dx, dx, Simplify]
Out[28]= (a3 dx^5)/6 + a0 x0^2 + 1/2 dx x0 (2 a0 + a1 x0) +
1/5 dx^4 (a2 + 2 a3 x0) + 1/3 dx^2 (a0 + x0 (2 a1 + a2 x0)) +
1/4 dx^3 (a1 + x0 (2 a2 + a3 x0))
In[29]:= r2 = Collect[(r1 - a0*x0^2)/dx, dx, Simplify]
Out[29]= (a3 dx^4)/6 + 1/2 x0 (2 a0 + a1 x0) +
1/5 dx^3 (a2 + 2 a3 x0) + 1/3 dx (a0 + x0 (2 a1 + a2 x0)) +
1/4 dx^2 (a1 + x0 (2 a2 + a3 x0))
In[30]:= r3 = Collect[(r2 - 1/2 x0 (2 a0 + a1 x0))/dx, dx, Simplify]
Out[30]= (a3 dx^3)/6 + 1/5 dx^2 (a2 + 2 a3 x0) +
1/3 (a0 + x0 (2 a1 + a2 x0)) + 1/4 dx (a1 + x0 (2 a2 + a3 x0))
In[31]:= r4 =
Collect[(r3 - 1/3 (a0 + x0 (2 a1 + a2 x0)))/dx, dx, Simplify]
Out[31]= (a3 dx^2)/6 + 1/5 dx (a2 + 2 a3 x0) +
1/4 (a1 + x0 (2 a2 + a3 x0))
In[32]:= r5 =
Collect[(r4 - 1/4 (a1 + x0 (2 a2 + a3 x0)))/dx, dx, Simplify]
Out[32]= (a3 dx)/6 + 1/5 (a2 + 2 a3 x0)
In[33]:= r6 = Collect[(r5 - 1/5 (a2 + 2 a3 x0))/dx, dx, Simplify]
Out[33]= a3/6
Definition at line 422 of file spline.hpp.
◆ integrate() [2/2]
|
inline |
Integrate spline with r^m prefactor.
Definition at line 579 of file spline.hpp.
◆ scale()
|
inline |
Definition at line 585 of file spline.hpp.
◆ coeffs() [1/2]
|
inline |
Definition at line 595 of file spline.hpp.
◆ coeffs() [2/2]
|
inline |
Definition at line 600 of file spline.hpp.
◆ values()
|
inline |
Definition at line 612 of file spline.hpp.
Member Data Documentation
◆ coeffs_
|
private |
Array of spline coefficients.
Definition at line 79 of file spline.hpp.
The documentation for this class was generated from the following file:
