Back to index

nux  3.0.0
Spline.h
Go to the documentation of this file.
00001 /*
00002  * Copyright 2010 Inalogic® Inc.
00003  *
00004  * This program is free software: you can redistribute it and/or modify it
00005  * under the terms of the GNU Lesser General Public License, as
00006  * published by the  Free Software Foundation; either version 2.1 or 3.0
00007  * of the License.
00008  *
00009  * This program is distributed in the hope that it will be useful, but
00010  * WITHOUT ANY WARRANTY; without even the implied warranties of
00011  * MERCHANTABILITY, SATISFACTORY QUALITY or FITNESS FOR A PARTICULAR
00012  * PURPOSE.  See the applicable version of the GNU Lesser General Public
00013  * License for more details.
00014  *
00015  * You should have received a copy of both the GNU Lesser General Public
00016  * License along with this program. If not, see <http://www.gnu.org/licenses/>
00017  *
00018  * Authored by: Jay Taoko <jaytaoko@inalogic.com>
00019  *
00020  */
00021 
00022 
00023 #ifndef SPLINE_H
00024 #define SPLINE_H
00025 
00026 #include <vector>
00027 
00028 namespace nux
00029 {
00030 
00031   double *d3_np_fs ( int n, double a[], double b[] );
00032   double *spline_cubic_set ( int n, double t[], double y[], int ibcbeg, double ybcbeg, int ibcend, double ybcend );
00033   double spline_cubic_val ( int n, double t[], double tval, double y[], double ypp[], double *ypval, double *yppval );
00034 
00035 
00036   class CubicSpline
00037   {
00038   public:
00039     CubicSpline();
00040     CubicSpline (const CubicSpline &Other);
00041     CubicSpline &operator = (const CubicSpline &Other);
00042 
00043     CubicSpline (int numpoint, std::vector<double> x_array, std::vector<double> y_array, int ibcbeg = 0, double ybcbeg = 0, int ibcend = 0, double ybcend = 0);
00044     void Set (int numpoint, std::vector<double> x_array, std::vector<double> y_array, int ibcbeg = 0, double ybcbeg = 0, int ibcend = 0, double ybcend = 0);
00045 
00046     ~CubicSpline();
00047 
00048     //**********************************************************************
00049     //  double* CubicSpline::Compute ( int n, double t[], double y[], int ibcbeg, double ybcbeg, int ibcend, double ybcend )
00050     //  Purpose:
00051     //
00052     //    SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline.
00053     //
00054     //  Discussion:
00055     //
00056     //    For data interpolation, the user must call SPLINE_SET to determine
00057     //    the second derivative data, passing in the data to be interpolated,
00058     //    and the desired boundary conditions.
00059     //
00060     //    The data to be interpolated, plus the SPLINE_SET output, defines
00061     //    the spline.  The user may then call SPLINE_VAL to evaluate the
00062     //    spline at any point.
00063     //
00064     //    The cubic spline is a piecewise cubic polynomial.  The intervals
00065     //    are determined by the "knots" or abscissas of the data to be
00066     //    interpolated.  The cubic spline has continous first and second
00067     //    derivatives over the entire interval of interpolation.
00068     //
00069     //    For any point T in the interval T(IVAL), T(IVAL+1), the form of
00070     //    the spline is
00071     //
00072     //      SPL(T) = A(IVAL)
00073     //             + B(IVAL) * ( T - T(IVAL) )
00074     //             + C(IVAL) * ( T - T(IVAL) )**2
00075     //             + D(IVAL) * ( T - T(IVAL) )**3
00076     //
00077     //    If we assume that we know the values Y(*) and DDY(*), which represent
00078     //    the values and second derivatives of the spline at each knot, then
00079     //    the coefficients can be computed as:
00080     //
00081     //      A(IVAL) = Y(IVAL)
00082     //      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
00083     //        - ( DDY(IVAL+1) + 2 * DDY(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
00084     //      C(IVAL) = DDY(IVAL) / 2
00085     //      D(IVAL) = ( DDY(IVAL+1) - DDY(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
00086     //
00087     //    Since the first derivative of the spline is
00088     //
00089     //      SPL'(T) =     B(IVAL)
00090     //              + 2 * C(IVAL) * ( T - T(IVAL) )
00091     //              + 3 * D(IVAL) * ( T - T(IVAL) )**2,
00092     //
00093     //    the requirement that the first derivative be continuous at interior
00094     //    knot I results in a total of N-2 equations, of the form:
00095     //
00096     //      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1)) + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))^2 = B(IVAL)
00097     //
00098     //    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
00099     //
00100     //      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) - ( DDY(IVAL) + 2 * DDY(IVAL-1) ) * H(IVAL-1) / 6 + DDY(IVAL-1) * H(IVAL-1)
00101     //      + ( DDY(IVAL) - DDY(IVAL-1) ) * H(IVAL-1) / 2
00102     //      = ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) - ( DDY(IVAL+1) + 2 * DDY(IVAL) ) * H(IVAL) / 6
00103     //
00104     //    or
00105     //
00106     //      DDY(IVAL-1) * H(IVAL-1) + 2 * DDY(IVAL) * ( H(IVAL-1) + H(IVAL) ) + DDY(IVAL) * H(IVAL)
00107     //      = 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
00108     //
00109     //    Boundary conditions must be applied at the first and last knots.
00110     //    The resulting tridiagonal system can be solved for the DDY values.
00111     //
00112     //
00113     //  Parameters:
00114     //
00115     //    Input, int N, the number of data points.  N must be at least 2.
00116     //    In the special case where N = 2 and IBCBEG = IBCEND = 0, the
00117     //    spline will actually be linear.
00118     //
00119     //    Input, double T[N], the knot values, that is, the points were data is
00120     //    specified.  The knot values should be distinct, and increasing.
00121     //
00122     //    Input, double Y[N], the data values to be interpolated.
00123     //
00124     //    Input, int IBCBEG, left boundary condition flag:
00125     //      0: the cubic spline should be a quadratic over the first interval;
00126     //      1: the first derivative at the left endpoint should be YBCBEG;
00127     //      2: the second derivative at the left endpoint should be YBCBEG.
00128     //
00129     //    Input, double YBCBEG, the values to be used in the boundary
00130     //    conditions if IBCBEG is equal to 1 or 2.
00131     //
00132     //    Input, int IBCEND, right boundary condition flag:
00133     //      0: the cubic spline should be a quadratic over the last interval;
00134     //      1: the first derivative at the right endpoint should be YBCEND;
00135     //      2: the second derivative at the right endpoint should be YBCEND.
00136     //
00137     //    Input, double YBCEND, the values to be used in the boundary
00138     //    conditions if IBCEND is equal to 1 or 2.
00139     //
00140     //    Output, double SPLINE_CUBIC_SET[N], the second derivatives of the cubic spline.
00141     //
00142     double *Compute (int ibcbeg, double ybcbeg, int ibcend, double ybcend);
00143 
00144     //**********************************************************************
00145     //  double Eval ( int n, double t[], double tval, double y[], double ddy[], double *dyval, double *ddyval )
00146     //  Purpose:
00147     //
00148     //    SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
00149     //
00150     //  Discussion:
00151     //
00152     //    SPLINE_CUBIC_SET must have already been called to define the values of YPP.
00153     //
00154     //    For any point T in the interval T(IVAL), T(IVAL+1), the form of
00155     //    the spline is
00156     //
00157     //      SPL(T) = A
00158     //             + B * ( T - T(IVAL) )
00159     //             + C * ( T - T(IVAL) )**2
00160     //             + D * ( T - T(IVAL) )**3
00161     //
00162     //    Here:
00163     //      A = Y(IVAL)
00164     //      B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
00165     //        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
00166     //      C = YPP(IVAL) / 2
00167     //      D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
00168     //
00169     //  Modified:
00170     //
00171     //    04 February 1999
00172     //
00173     //  Author:
00174     //
00175     //    John Burkardt
00176     //
00177     //  Parameters:
00178     //
00179     //    Input, int N, the number of knots.
00180     //
00181     //    Input, double Y[N], the data values at the knots.
00182     //
00183     //    Input, double T[N], the knot values.
00184     //
00185     //    Input, double TVAL, a point, typically between T[0] and T[N-1], at
00186     //    which the spline is to be evalulated.  If TVAL lies outside
00187     //    this range, extrapolation is used.
00188     //
00189     //    Input, double Y[N], the data values at the knots.
00190     //
00191     //    Input, double YPP[N], the second derivatives of the spline at
00192     //    the knots.
00193     //
00194     //    Output, double *YPVAL, the derivative of the spline at TVAL.
00195     //
00196     //    Output, double *YPPVAL, the second derivative of the spline at TVAL.
00197     //
00198     //    Output, double SPLINE_VAL, the value of the spline at TVAL.
00199     //
00200     double Eval (double tval);
00201 
00202     //**********************************************************************
00203     //  double* SolveTridiag ( int n, double a[], double b[] );
00204     //  Purpose:
00205     //
00206     //    D3_NP_FS factors and solves a D3 system.
00207     //
00208     //  Discussion:
00209     //
00210     //    The D3 storage format is used for a tridiagonal matrix.
00211     //    The superdiagonal is stored in entries (1,2:N), the diagonal in
00212     //    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
00213     //    original matrix is "collapsed" vertically into the array.
00214     //
00215     //    This algorithm requires that each diagonal entry be nonzero.
00216     //    It does not use pivoting, and so can fail on systems that
00217     //    are actually nonsingular.
00218     //
00219     //  Example:
00220     //
00221     //    Here is how a D3 matrix of order 5 would be stored:
00222     //
00223     //       *  A12 A23 A34 A45
00224     //      A11 A22 A33 A44 A55
00225     //      A21 A32 A43 A54  *
00226     //
00227     //  Parameters:
00228     //
00229     //    Input, int N, the order of the linear system.
00230     //
00231     //    Input/output, double A[3*N].
00232     //    On input, the nonzero diagonals of the linear system.
00233     //    On output, the data in these vectors has been overwritten
00234     //    by factorization information.
00235     //
00236     //    Input, double B[N], the right hand side.
00237     //
00238     //    Output, double D3_NP_FS[N], the solution of the linear system.
00239     //    This is NULL if there was an error because one of the diagonal
00240     //    entries was zero.
00241     //
00242     double *SolveTridiag ( int n, double a[], double b[] );
00243 
00244   public:
00245     double *t;
00246     double *y;
00247     double *ddy;
00248 
00249     int ibcbeg_;
00250     double ybcbeg_;
00251     int ibcend_;
00252     double ybcend_;
00253 
00254     int np;  // number of points
00255   };
00256 
00257 }
00258 
00259 #endif // SPLINE_H