1/*! 2 * \file include/TFEL/Math/CubicSpline.ixx 3 * \author Castelier Etienne 4 * \date 07/08/2007 5 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights 6 * reserved. 7 * This project is publicly released under either the GNU GPL Licence 8 * or the CECILL-A licence. A copy of thoses licences are delivered 9 * with the sources of TFEL. CEA or EDF may also distribute this 10 * project under specific licensing conditions. 11 */ 12 13#ifndef LIB_TFEL_MATH_CUBICSPLINE_IXX 14#define LIB_TFEL_MATH_CUBICSPLINE_IXX 1 15 16#include <limits> 17#include <algorithm> 18 19#include "TFEL/Math/General/Abs.hxx" 20 21namespace tfel { 22 23 namespace math { 24 25 template <typename real, typename value> 26 bool CubicSpline<real, value>::PointComparator::operator()( 27 const typename CubicSpline<real, value>::Point& p, 28 const real& x) const { 29 return p.x < x; 30 } // end of CubicSpline<real,value>::Point 31 32 template <typename real, typename value> 33 template <typename AIterator, typename OIterator> 34 void CubicSpline<real, value>::setCollocationPoints(AIterator px, 35 AIterator pxe, 36 OIterator py) { 37 using namespace std; 38 AIterator px2; 39 if (px == pxe) { 40 throw(CubicSplineInvalidAbscissaVectorSize()); 41 } 42 this->values.clear(); 43 Point p; 44 p.x = *px; 45 p.y = *py; 46 this->values.push_back(p); 47 px2 = px; 48 ++px; 49 ++py; 50 while (px != pxe) { 51 px2 = px - 1u; 52 if (*px2 >= *px) { 53 throw(CubicSplineUnorderedAbscissaVector()); 54 } 55 p.x = *px; 56 p.y = *py; 57 this->values.push_back(p); 58 px2 = px; 59 ++px; 60 ++py; 61 } 62 this->buildInterpolation(); 63 } // CubicSpline<real,value>::CubicSpline 64 65 template <typename real, typename value> 66 template <typename AContainer, typename OContainer> 67 void CubicSpline<real, value>::setCollocationPoints(const AContainer& x, 68 const OContainer& y) { 69 if (x.size() < 1) { 70 throw(CubicSplineInvalidAbscissaVectorSize()); 71 } 72 if (y.size() < 1) { 73 throw(CubicSplineInvalidOrdinateVectorSize()); 74 } 75 if (x.size() != y.size()) { 76 throw(CubicSplineInvalidInputs()); 77 } 78 this->setCollocationPoints(x.begin(), x.end(), y.begin()); 79 } // CubicSpline<real,value>::CubicSpline 80 81 template <typename real, typename value> 82 void CubicSpline<real, value>::buildInterpolation() { 83 using namespace std; 84 using std::vector; 85 if (this->values.size() == 1) { 86 this->values[0].d = value(real(0)); 87 } else { 88 typename vector<real>::size_type i; 89 typename vector<real>::size_type s = this->values.size() - 1u; 90 vector<real> m(2 * s + 1u); // matrix (main an upper diagonal) 91 real* const md = &m[0]; // main diagonal 92 real* const mu = md + s + 1u; // upper diagonal 93 real ho = real(0); 94 value uo = value(real(0)); 95 for (i = 0; i != s; ++i) { 96 const real hn = 1. / (this->values[i + 1].x - this->values[i].x); 97 const value un = 98 3. * hn * hn * (this->values[i + 1].y - this->values[i].y); 99 mu[i] = hn; 100 md[i] = 2. * (hn + ho); 101 this->values[i].d = un + uo; 102 uo = un; 103 ho = hn; 104 } 105 md[s] = 2. * ho; 106 this->values[s].d = uo; 107 solveTridiagonalLinearSystem(mu, md); 108 } 109 } 110 111 template <typename real, typename value> 112 void CubicSpline<real, value>::solveTridiagonalLinearSystem( 113 const real* const c, real* const b) { 114 using namespace std; 115 using std::vector; 116 typename vector<typename CubicSpline<real, value>::Point>::size_type i; 117 typename vector<typename CubicSpline<real, value>::Point>::size_type n; 118 const real prec = 100 * numeric_limits<real>::min(); 119 n = this->values.size(); 120 for (i = 1; i < n; i++) { 121 if (abs(b[i - 1]) < prec) { 122 throw(CubicSplineNullPivot()); 123 } 124 real m = c[i - 1] / b[i - 1]; 125 b[i] -= m * c[i - 1]; 126 this->values[i].d -= m * this->values[i - 1].d; 127 } 128 if (abs(b[n - 1]) < prec) { 129 throw(CubicSplineNullPivot()); 130 } 131 this->values[n - 1].d /= b[n - 1]; 132 i = n; 133 i -= 2u; 134 for (; i != 0; i--) { 135 if (abs(b[i]) < prec) { 136 throw(CubicSplineNullPivot()); 137 } 138 this->values[i].d = 139 (this->values[i].d - c[i] * (this->values[i + 1].d)) / b[i]; 140 } 141 if (abs(b[0]) < prec) { 142 throw(CubicSplineNullPivot()); 143 } 144 this->values[0].d = 145 (this->values[0].d - c[0] * (this->values[1].d)) / b[0]; 146 } 147 148 template <typename real, typename value> 149 void CubicSpline<real, value>::computeLocalCoefficients( 150 value& a2, 151 value& a3, 152 const typename std::vector< 153 typename CubicSpline<real, value>::Point>::const_iterator pa) { 154 using namespace std; 155 const typename vector<Point>::const_iterator pb = pa + 1; 156 const real usL = 1 / (pb->x - pa->x); 157 const value Dy = (pb->y - pa->y) * usL; 158 a2 = (3 * Dy - pb->d - 2 * pa->d) * usL; 159 a3 = (-2 * Dy + pb->d + pa->d) * usL * usL; 160 } // end of CubicSpline<real,value>::computeLocalCoefficient 161 162 template <typename real, typename value> 163 value CubicSpline<real, value>::computeLocalIntegral( 164 const real xa, 165 const real xb, 166 const typename std::vector< 167 typename CubicSpline<real, value>::Point>::const_iterator pa) { 168 using namespace std; 169 const value py = pa->y; 170 const value d = pa->d; 171 value a2; 172 value a3; 173 CubicSpline<real, value>::computeLocalCoefficients(a2, a3, pa); 174 const real x0 = xa - pa->x; 175 const real x1 = xb - pa->x; 176 return (3 * a3 * (x1 * x1 * x1 * x1 - x0 * x0 * x0 * x0) + 177 4 * a2 * (x1 * x1 * x1 - x0 * x0 * x0) + 178 6 * d * (x1 * x1 - x0 * x0) + 12 * py * (x1 - x0)) / 179 12; 180 } // end of 181 182 template <typename real, typename value> 183 value CubicSpline<real, value>::computeMeanValue(const real xa, 184 const real xb) const { 185 return this->computeIntegral(xa, xb) / (xb - xa); 186 } // end of value CubicSpline<real,value>::computeMeanValue 187 188 template <typename real, typename value> 189 value CubicSpline<real, value>::computeIntegral(const real xa, 190 const real xb) const { 191 using namespace std; 192 typename vector<Point>::const_iterator pa; 193 typename vector<Point>::const_iterator pb; 194 if (this->values.empty()) { 195 throw(CubicSplineUninitialised()); 196 } 197 if (this->values.size() == 1) { 198 const auto& f = values.back().y; 199 return f * (xb - xa); 200 } 201 if (xb < xa) { 202 return -this->computeIntegral(xb, xa); 203 } 204 pa = lower_bound(this->values.begin(), this->values.end(), xa, 205 PointComparator()); 206 pb = lower_bound(this->values.begin(), this->values.end(), xb, 207 PointComparator()); 208 if (pa == pb) { 209 if (pb == this->values.begin()) { 210 const real xe = this->values.front().x; 211 const auto& ye = this->values.front().y; 212 const auto& df = this->values.front().d; 213 // l'équation de l'extrapoltion est : 214 // y = ye-df*xe + df*x 215 // l'intégrale est alors: 216 // (ye-df*xe)*(xb-xa)+0.5*df*(xb-xa)*(xb-xa) 217 return ye * (xb - xa) + 218 0.5 * df * ((xb - xe) * (xb - xe) - (xa - xe) * (xa - xe)); 219 } else if (pb == this->values.end()) { 220 const real xe = this->values.back().x; 221 const auto& ye = this->values.back().y; 222 const auto& df = this->values.back().d; 223 // l'équation de l'extrapoltion est : 224 // y = ye-df*xe + df*x 225 // l'intégrale est alors: 226 // (ye-df*xe)*(xb-xa)+0.5*df*(xb-xa)*(xb-xa) 227 return ye * (xb - xa) + 228 0.5 * df * ((xb - xe) * (xb - xe) - (xa - xe) * (xa - xe)); 229 } else { 230 pa = pb - 1; 231 // spline=pa->y+(x-pa->x)*(pa->d+(x-pa->x)*(a2+(x-pa->x)*a3)); 232 // =pa->y+(x-pa->x)*(pa->d+(x-pa->x)*(a2+(x-pa->x)*a3)); 233 return CubicSpline<real, value>::computeLocalIntegral(xa, xb, pa); 234 } 235 } 236 value s(real(0)); 237 if (pa == this->values.begin()) { 238 const real xe = this->values.front().x; 239 const auto& ye = this->values.front().y; 240 const auto& df = this->values.front().d; 241 // l'équation de l'extrapoltion est : 242 // y = ye-df*xe + df*x 243 // l'intégrale est alors: 244 // (ye-df*xe)*(pa->a-xa)+0.5*df*(pa->a-xa)*(pa->a-xa) 245 s += ye * (xe - xa) - 0.5 * df * ((xa - xe) * (xa - xe)); 246 } else { 247 s += CubicSpline<real, value>::computeLocalIntegral(xa, pa->x, pa - 1); 248 } 249 if (pb == this->values.end()) { 250 const real xe = this->values.back().x; 251 const auto& ye = this->values.back().y; 252 const auto& df = this->values.back().d; 253 // l'équation de l'extrapoltion est : 254 // y = ye-df*xe + df*x 255 // l'intégrale est alors: 256 // (ye-df*xe)*(xb-xa)+0.5*df*(xb-xa)*(xb-xa) 257 s += ye * (xb - xe) + 0.5 * df * ((xb - xe) * (xb - xe)); 258 // s += (ye-df*xe)*(xb-pb->x)+0.5*df*(xb-pb->x)*(xb-pb->x); 259 } else { 260 s += CubicSpline<real, value>::computeLocalIntegral((pb - 1)->x, xb, 261 pb - 1); 262 } 263 --pb; 264 while (pa != pb) { 265 s += CubicSpline<real, value>::computeLocalIntegral(pa->x, (pa + 1)->x, 266 pa); 267 ++pa; 268 } 269 return s; 270 } // end of CubicSpline<real,value>::computeIntegral 271 272 // Valeur de la spline et de sa dérivée 273 // Interpolation cubique 274 // x : Abscisse où calculer la valeur de la spline 275 // f : Valeur de la fonction 276 // df : Valeur de la dérivée 277 template <typename real, typename value> 278 void CubicSpline<real, value>::getValues(value& f, 279 value& df, 280 real x) const { 281 using namespace std; 282 typename vector<Point>::const_iterator in; 283 if (this->values.empty()) { 284 throw(CubicSplineUninitialised()); 285 } 286 // Cas d'un seul couple 287 if (this->values.size() == 1) { 288 f = this->values[0].y; 289 df = value(real(0)); 290 return; 291 } 292 in = lower_bound(this->values.begin(), this->values.end(), x, 293 PointComparator()); 294 // Extrapolation 295 if (in == this->values.begin()) { 296 x -= in->x; 297 df = in->d; 298 f = in->y + x * df; 299 return; 300 } 301 if (in == this->values.end()) { 302 --in; 303 x -= in->x; 304 df = in->d; 305 f = in->y + x * df; 306 return; 307 } 308 typename vector<Point>::const_iterator ip = in - 1; 309 value a2; 310 value a3; 311 CubicSpline<real, value>::computeLocalCoefficients(a2, a3, ip); 312 x -= ip->x; 313 f = ip->y + x * (ip->d + x * (a2 + x * a3)); 314 df = ip->d + x * (2 * a2 + x * 3 * a3); 315 } 316 317 template <typename real, typename value> 318 value CubicSpline<real, value>::operator()(const real x) const { 319 return this->getValue(x); 320 } 321 322 // Valeur de la spline 323 // Interpolation cubique 324 // x : Abscisse où calculer la valeur de la spline 325 template <typename real, typename value> 326 value CubicSpline<real, value>::getValue(real x) const { 327 using namespace std; 328 typename vector<Point>::const_iterator in; 329 if (this->values.empty()) { 330 throw(CubicSplineUninitialised()); 331 } 332 // Cas d'un seul couple 333 if (this->values.size() == 1) { 334 return this->values[0].y; 335 } 336 in = lower_bound(this->values.begin(), this->values.end(), x, 337 PointComparator()); 338 // Extrapolation 339 if (in == this->values.begin()) { 340 x -= in->x; 341 return in->y + x * in->d; 342 } 343 if (in == this->values.end()) { 344 --in; 345 x -= in->x; 346 return in->y + x * in->d; 347 } 348 typename vector<Point>::const_iterator ip = in - 1; 349 value a2; 350 value a3; 351 CubicSpline<real, value>::computeLocalCoefficients(a2, a3, ip); 352 x -= ip->x; 353 return ip->y + x * (ip->d + x * (a2 + x * a3)); 354 } 355 356 // Valeur de la spline et de sa dérivée 357 // Interpolation cubique 358 // x : Abscisse où calculer la valeur de la spline 359 // f : Valeur de la fonction 360 // df : Valeur de la dérivée 361 // d2f : Valeur de la dérivée seconde 362 template <typename real, typename value> 363 void CubicSpline<real, value>::getValues(value& f, 364 value& df, 365 value& d2f, 366 real x) const { 367 using namespace std; 368 if (this->values.empty()) { 369 throw(CubicSplineUninitialised()); 370 } 371 typename vector<Point>::const_iterator in; 372 in = this->values.begin(); 373 // Cas d'un seul couple 374 if (this->values.size() == 1) { 375 f = this->values[0].y; 376 df = value(real(0)); 377 d2f = value(real(0)); 378 return; 379 } 380 in = lower_bound(this->values.begin(), this->values.end(), x, 381 PointComparator()); 382 // Extrapolation 383 if (in == this->values.begin()) { 384 x -= in->x; 385 df = in->d; 386 f = in->y + x * df; 387 d2f = value(real(0)); 388 return; 389 } 390 if (in == this->values.end()) { 391 --in; 392 x -= in->x; 393 df = in->d; 394 f = in->y + x * df; 395 d2f = value(real(0)); 396 return; 397 } 398 typename vector<Point>::const_iterator ip = in - 1; 399 value a2; 400 value a3; 401 CubicSpline<real, value>::computeLocalCoefficients(a2, a3, ip); 402 x -= ip->x; 403 f = ip->y + x * (ip->d + x * (a2 + x * a3)); 404 df = ip->d + x * (2 * a2 + x * 3 * a3); 405 d2f = 2 * a2 + x * 6 * a3; 406 } 407 408 } // end of namespace math 409 410} // end of namespace tfel 411 412#endif /* LIB_TFEL_MATH_CUBICSPLINE_IXX */ 413