nux-1.16.0
|
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