ADMB Documentation  11.1.2274
 All Classes Files Functions Variables Typedefs Friends Defines
cspline.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: cspline.cpp 1711 2014-02-28 22:46:44Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00012 #include <fvar.hpp>
00013 
00014 dvector spline(const dvector &x, const dvector&y, double yp1,double ypn);
00015 dvector spline_cubic_set(int n, const dvector& t, const dvector& y, int ibcbeg,
00016                          double ybcbeg, int ibcend, double ybcend );
00017 double spline_cubic_val(int n, const dvector& t, double tval,
00018                         const dvector& y, const dvector& ypp);
00019 
00020 double splint(const dvector& xa, const dvector& ya, const dvector& y2a,
00021   double x);
00022 
00030 cubic_spline_function::cubic_spline_function(const dvector & _x,
00031   const dvector& _y, double yp1, double ypn) : x(_x) , y(_y)
00032 {
00033   y2.allocate(x);
00034   y2=spline(x,y,yp1,ypn);
00035 }
00036 
00041 double cubic_spline_function::operator () (double u)
00042 {
00043   // need to deal with u<x(1) or u>x(2)
00044   return splint(x,y,y2,u);
00045 }
00046 
00051 dvector cubic_spline_function::operator()(const dvector& u)
00052 {
00053   int mmin=u.indexmin();
00054   int mmax=u.indexmax();
00055   dvector z(mmin,mmax);
00056   for (int i=mmin;i<=mmax;i++)
00057   {
00058     z(i)=splint(x,y,y2,u(i));
00059   }
00060   return z;
00061 }
00062 
00074 dvector spline(const dvector& _x, const dvector& _y,double yp1,double ypn)
00075 {
00076   int ibcbeg, ibcend;
00077   double ybcbeg, ybcend;
00078   dvector x = _x;
00079   x.shift(0);
00080   dvector y = _y;
00081   y.shift(0);
00082   if(yp1 > 0.99e30 )
00083   {
00084     ibcbeg = 2;
00085     ybcbeg = 0.0;
00086   }
00087   else
00088   {
00089     ibcbeg = 1;
00090     ybcbeg = yp1;
00091   }
00092   if(ypn > 0.99e30 )
00093   {
00094     ibcend = 2;
00095     ybcend = 0.0;
00096   }
00097   else
00098   {
00099     ibcend = 1;
00100     ybcend = ypn;
00101   }
00102   dvector ret = spline_cubic_set(x.size(), x, y, ibcbeg, ybcbeg, ibcend,
00103     ybcend);
00104   ret.shift(_x.indexmin());
00105   return ret;
00106 }
00107 
00117 double splint(const dvector& _xa, const dvector& _ya, const dvector& _y2a,
00118   double x)
00119 {
00120   return spline_cubic_val(_xa.size(), _xa, x, _ya, _y2a);
00121 }
00122 
00135 double spline_cubic_val(int n, const dvector& t, double tval,
00136                         const dvector& y, const dvector& ypp)
00137 //
00138 //  Purpose:
00139 //
00140 //    SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
00141 //
00142 //  Discussion:
00143 //
00144 //    SPLINE_CUBIC_SET must have already been called to define the values of
00145 //      YPP.
00146 //
00147 //    For any point T in the interval T(IVAL), T(IVAL+1), the form of
00148 //    the spline is
00149 //
00150 //      SPL(T) = A
00151 //             + B * ( T - T(IVAL) )
00152 //             + C * ( T - T(IVAL) )**2
00153 //             + D * ( T - T(IVAL) )**3
00154 //
00155 //    Here:
00156 //      A = Y(IVAL)
00157 //      B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
00158 //        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
00159 //      C = YPP(IVAL) / 2
00160 //      D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
00161 //
00162 //  Licensing:
00163 //
00164 //    This code is distributed under the GNU LGPL license.
00165 //
00166 //  Modified:
00167 //
00168 //    04 February 1999
00169 //
00170 //  Author:
00171 //
00172 //    John Burkardt
00173 //
00174 //  Modified By:
00175 //
00176 //    Derek Seiple (20 August 2010)
00177 //
00178 //  Parameters:
00179 //
00180 //    Input, int N, the number of knots.
00181 //
00182 //    Input, double Y[N], the data values at the knots.
00183 //
00184 //    Input, double T[N], the knot values.
00185 //
00186 //    Input, double TVAL, a point, typically between T[0] and T[N-1], at
00187 //    which the spline is to be evalulated.  If TVAL lies outside
00188 //    this range, extrapolation is used.
00189 //
00190 //    Input, double YPP[N], the second derivatives of the spline at
00191 //    the knots.
00192 //
00193 //    Output, double SPLINE_VAL, the value of the spline at TVAL.
00194 //
00195 {
00196   double dt;
00197   double h;
00198   int i;
00199   int ival;
00200   double yval;
00201 
00202   //  Determine the interval [ T(I), T(I+1) ] that contains TVAL.
00203   //  Values below T[0] or above T[N-1] use extrapolation.
00204   ival = n - 2;
00205 
00206   for ( i = 0; i < n-1; i++ )
00207   {
00208     if ( tval < t[i+1] )
00209     {
00210       ival = i;
00211       break;
00212     }
00213   }
00214 
00215   //  In the interval I, the polynomial is in terms of a normalized
00216   //  coordinate between 0 and 1.
00217   dt = tval - t[ival];
00218   h = t[ival+1] - t[ival];
00219 
00220   yval = y[ival]
00221     + dt * ( ( y[ival+1] - y[ival] ) / h
00222            - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
00223     + dt * ( 0.5 * ypp[ival]
00224     + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) );
00225 
00226   /* *ypval = ( y[ival+1] - y[ival] ) / h
00227     - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
00228     + dt * ( ypp[ival]
00229     + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) );
00230 
00231   *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h;*/
00232 
00233   return yval;
00234 }
00235 
00244 double *d3_np_fs(int n, const dvector& _a, const dvector& _b)
00245 //
00246 //  Purpose:
00247 //
00248 //    D3_NP_FS factors and solves a D3 system.
00249 //
00250 //  Discussion:
00251 //
00252 //    The D3 storage format is used for a tridiagonal matrix.
00253 //    The superdiagonal is stored in entries (1,2:N), the diagonal in
00254 //    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
00255 //    original matrix is "collapsed" vertically into the array.
00256 //
00257 //    This algorithm requires that each diagonal entry be nonzero.
00258 //    It does not use pivoting, and so can fail on systems that
00259 //    are actually nonsingular.
00260 //
00261 //  Example:
00262 //
00263 //    Here is how a D3 matrix of order 5 would be stored:
00264 //
00265 //       *  A12 A23 A34 A45
00266 //      A11 A22 A33 A44 A55
00267 //      A21 A32 A43 A54  *
00268 //
00269 //  Licensing:
00270 //
00271 //    This code is distributed under the GNU LGPL license.
00272 //
00273 //  Modified:
00274 //
00275 //    15 November 2003
00276 //
00277 //  Author:
00278 //
00279 //    John Burkardt
00280 //
00281 //  Modified By:
00282 //
00283 //    Derek Seiple (20 August 2010)
00284 //
00285 //  Parameters:
00286 //
00287 //    Input, int N, the order of the linear system.
00288 //
00289 //    Input/output, double A[3*N].
00290 //    On input, the nonzero diagonals of the linear system.
00291 //    On output, the data in these vectors has been overwritten
00292 //    by factorization information.
00293 //
00294 //    Input, double B[N], the right hand side.
00295 //
00296 //    Output, double D3_NP_FS[N], the solution of the linear system.
00297 //    This is NULL if there was an error because one of the diagonal
00298 //    entries was zero.
00299 //
00300 {
00301   ADUNCONST(dvector,a)
00302   ADUNCONST(dvector,b)
00303   int i;
00304   double *x;
00305   double xmult;
00306 
00307   //  Check.
00308   for ( i = 0; i < n; i++ )
00309   {
00310     if ( a[1+i*3] == 0.0 )
00311     {
00312       return NULL;
00313     }
00314   }
00315   x = new double[n];
00316 
00317   for ( i = 0; i < n; i++ )
00318   {
00319     x[i] = b[i];
00320   }
00321 
00322   for ( i = 1; i < n; i++ )
00323   {
00324     xmult = a[2+(i-1)*3] / a[1+(i-1)*3];
00325     a[1+i*3] = a[1+i*3] - xmult * a[0+i*3];
00326     x[i] = x[i] - xmult * x[i-1];
00327   }
00328 
00329   x[n-1] = x[n-1] / a[1+(n-1)*3];
00330   for ( i = n-2; 0 <= i; i-- )
00331   {
00332     x[i] = ( x[i] - a[0+(i+1)*3] * x[i+1] ) / a[1+i*3];
00333   }
00334 
00335   return x;
00336 }
00337 
00358 dvector spline_cubic_set(int n, const dvector& t, const dvector& y, int ibcbeg,
00359   double ybcbeg, int ibcend, double ybcend )
00360 //
00361 //  Purpose:
00362 //
00363 //    SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic
00364 //      spline.
00365 //
00366 //  Discussion:
00367 //
00368 //    For data interpolation, the user must call SPLINE_SET to determine
00369 //    the second derivative data, passing in the data to be interpolated,
00370 //    and the desired boundary conditions.
00371 //
00372 //    The data to be interpolated, plus the SPLINE_SET output, defines
00373 //    the spline.  The user may then call SPLINE_VAL to evaluate the
00374 //    spline at any point.
00375 //
00376 //    The cubic spline is a piecewise cubic polynomial.  The intervals
00377 //    are determined by the "knots" or abscissas of the data to be
00378 //    interpolated.  The cubic spline has continous first and second
00379 //    derivatives over the entire interval of interpolation.
00380 //
00381 //    For any point T in the interval T(IVAL), T(IVAL+1), the form of
00382 //    the spline is
00383 //
00384 //      SPL(T) = A(IVAL)
00385 //             + B(IVAL) * ( T - T(IVAL) )
00386 //             + C(IVAL) * ( T - T(IVAL) )**2
00387 //             + D(IVAL) * ( T - T(IVAL) )**3
00388 //
00389 //    If we assume that we know the values Y(*) and YPP(*), which represent
00390 //    the values and second derivatives of the spline at each knot, then
00391 //    the coefficients can be computed as:
00392 //
00393 //      A(IVAL) = Y(IVAL)
00394 //      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
00395 //        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
00396 //      C(IVAL) = YPP(IVAL) / 2
00397 //      D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
00398 //
00399 //    Since the first derivative of the spline is
00400 //
00401 //      SPL'(T) =     B(IVAL)
00402 //              + 2 * C(IVAL) * ( T - T(IVAL) )
00403 //              + 3 * D(IVAL) * ( T - T(IVAL) )**2,
00404 //
00405 //    the requirement that the first derivative be continuous at interior
00406 //    knot I results in a total of N-2 equations, of the form:
00407 //
00408 //      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
00409 //      + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)
00410 //
00411 //    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
00412 //
00413 //      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
00414 //      - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
00415 //      + YPP(IVAL-1) * H(IVAL-1)
00416 //      + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
00417 //      =
00418 //      ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
00419 //      - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
00420 //
00421 //    or
00422 //
00423 //      YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
00424 //      + YPP(IVAL) * H(IVAL)
00425 //      =
00426 //      6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
00427 //      - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
00428 //
00429 //    Boundary conditions must be applied at the first and last knots.
00430 //    The resulting tridiagonal system can be solved for the YPP values.
00431 //
00432 //  Licensing:
00433 //
00434 //    This code is distributed under the GNU LGPL license.
00435 //
00436 //  Modified:
00437 //
00438 //    06 February 2004
00439 //
00440 //  Author:
00441 //
00442 //    John Burkardt
00443 //
00444 //  Modified By:
00445 //
00446 //    Derek Seiple (20 August 2010)
00447 //
00448 //  Parameters:
00449 //
00450 //    Input, int N, the number of data points.  N must be at least 2.
00451 //    In the special case where N = 2 and IBCBEG = IBCEND = 0, the
00452 //    spline will actually be linear.
00453 //
00454 //    Input, double T[N], the knot values, that is, the points were data is
00455 //    specified.  The knot values should be distinct, and increasing.
00456 //
00457 //    Input, double Y[N], the data values to be interpolated.
00458 //
00459 //    Input, int IBCBEG, left boundary condition flag:
00460 //      0: the cubic spline should be a quadratic over the first interval;
00461 //      1: the first derivative at the left endpoint should be YBCBEG;
00462 //      2: the second derivative at the left endpoint should be YBCBEG.
00463 //
00464 //    Input, double YBCBEG, the values to be used in the boundary
00465 //    conditions if IBCBEG is equal to 1 or 2.
00466 //
00467 //    Input, int IBCEND, right boundary condition flag:
00468 //      0: the cubic spline should be a quadratic over the last interval;
00469 //      1: the first derivative at the right endpoint should be YBCEND;
00470 //      2: the second derivative at the right endpoint should be YBCEND.
00471 //
00472 //    Input, double YBCEND, the values to be used in the boundary
00473 //    conditions if IBCEND is equal to 1 or 2.
00474 //
00475 //    Output, double SPLINE_CUBIC_SET[N], the second derivatives of the cubic
00476 //      spline.
00477 //
00478 {
00479   dvector a(0,3*n-1);
00480   dvector b(0,n-1);
00481   int i;
00482   double *ypp;
00483 
00484   //  Check.
00485   if ( n <= 1 )
00486   {
00487     cout << "\n";
00488     cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00489     cout << "  The number of data points N must be at least 2.\n";
00490     cout << "  The input value is " << n << ".\n";
00491     //return NULL;
00492   }
00493 
00494   for ( i = 0; i < n - 1; i++ )
00495   {
00496     if ( t[i+1] <= t[i] )
00497     {
00498       cout << "\n";
00499       cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00500       cout << "  The knots must be strictly increasing, but\n";
00501       cout << "  T(" << i   << ") = " << t[i]   << "\n";
00502       cout << "  T(" << i+1 << ") = " << t[i+1] << "\n";
00503     }
00504   }
00505 
00506   //  Set up the first equation.
00507   if ( ibcbeg == 0 )
00508   {
00509     b[0] = 0.0;
00510     a[1+0*3] = 1.0;
00511     a[0+1*3] = -1.0;
00512   }
00513   else if ( ibcbeg == 1 )
00514   {
00515     b[0] = ( y[1] - y[0] ) / ( t[1] - t[0] ) - ybcbeg;
00516     a[1+0*3] = ( t[1] - t[0] ) / 3.0;
00517     a[0+1*3] = ( t[1] - t[0] ) / 6.0;
00518   }
00519   else if ( ibcbeg == 2 )
00520   {
00521     b[0] = ybcbeg;
00522     a[1+0*3] = 1.0;
00523     a[0+1*3] = 0.0;
00524   }
00525   else
00526   {
00527     cout << "\n";
00528     cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00529     cout << "  IBCBEG must be 0, 1 or 2.\n";
00530     cout << "  The input value is " << ibcbeg << ".\n";
00531   }
00532 
00533   //  Set up the intermediate equations.
00534   for ( i = 1; i < n-1; i++ )
00535   {
00536     b[i] = ( y[i+1] - y[i] ) / ( t[i+1] - t[i] )
00537       - ( y[i] - y[i-1] ) / ( t[i] - t[i-1] );
00538     a[2+(i-1)*3] = ( t[i] - t[i-1] ) / 6.0;
00539     a[1+ i   *3] = ( t[i+1] - t[i-1] ) / 3.0;
00540     a[0+(i+1)*3] = ( t[i+1] - t[i] ) / 6.0;
00541   }
00542 
00543   //  Set up the last equation.
00544   if ( ibcend == 0 )
00545   {
00546     b[n-1] = 0.0;
00547     a[2+(n-2)*3] = -1.0;
00548     a[1+(n-1)*3] = 1.0;
00549   }
00550   else if ( ibcend == 1 )
00551   {
00552     b[n-1] = ybcend - ( y[n-1] - y[n-2] ) / ( t[n-1] - t[n-2] );
00553     a[2+(n-2)*3] = ( t[n-1] - t[n-2] ) / 6.0;
00554     a[1+(n-1)*3] = ( t[n-1] - t[n-2] ) / 3.0;
00555   }
00556   else if ( ibcend == 2 )
00557   {
00558     b[n-1] = ybcend;
00559     a[2+(n-2)*3] = 0.0;
00560     a[1+(n-1)*3] = 1.0;
00561   }
00562   else
00563   {
00564     cout << "\n";
00565     cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00566     cout << "  IBCEND must be 0, 1 or 2.\n";
00567     cout << "  The input value is " << ibcend << ".\n";
00568   }
00569 
00570   //  Solve the linear system.
00571   if ( n == 2 && ibcbeg == 0 && ibcend == 0 )
00572   {
00573     dvector ret(0,1);
00574     ret[0] = 0.0;
00575     ret[1] = 0.0;
00576     return ret;
00577   }
00578   else
00579   {
00580     ypp = d3_np_fs ( n, a, b );
00581     dvector ret(0,n-1);
00582     if ( !ypp )
00583     {
00584       cout << "\n";
00585       cout << "SPLINE_CUBIC_SET - Fatal error!\n";
00586       cout << "  The linear system could not be solved.\n";
00587     }
00588     for(i=0;i<n;i++)
00589     {
00590        ret(i) =ypp[i];
00591     }
00592     return ret;
00593   }
00594 }