tools::spline::quintic Class Reference

Inheritance diagram for tools::spline::quintic:

Collaboration diagram for tools::spline::quintic:

Public Member Functions
	quintic (std::ostream &a_out, size_t a_n, double a_x[], double a_y[])
	quintic (const quintic &a_from)
quintic &	operator= (const quintic &a_from)
double	eval (double x) const
Public Member Functions inherited from tools::spline::base_spline
	base_spline (std::ostream &a_out, double delta, double xmin, double xmax, size_t np, bool step)
virtual	~base_spline ()

Protected Member Functions
	quintic (std::ostream &a_out)
size_t	find_x (double x) const
void	build_coeff ()
Protected Member Functions inherited from tools::spline::base_spline
	base_spline (std::ostream &a_out)
	base_spline (const base_spline &a_from)
base_spline &	operator= (const base_spline &a_from)

Protected Attributes
std::vector< quintic_poly >	fPoly
Protected Attributes inherited from tools::spline::base_spline
std::ostream &	m_out
double	fDelta
double	fXmin
double	fXmax
size_t	fNp
bool	fKstep

Detailed Description

Definition at line 416 of file spline.

Constructor & Destructor Documentation

quintic() [1/3]

tools::spline::quintic::quintic ( std::ostream &  a_out )

inlineprotected

Definition at line 418 of file spline.

:base_spline(a_out),fPoly() {}

quintic() [2/3]

tools::spline::quintic::quintic ( std::ostream &  a_out, size_t  a_n, double  a_x[], double  a_y[]  )

inline

Definition at line 420 of file spline.

  :base_spline(a_out,-1,0,0,a_n,false) {
    if(!a_n) {
      m_out << "tools::spline::quintic : a_np is null." << std::endl;
      return;
    }
    fXmin = a_x[0];
    fXmax = a_x[a_n-1];
    fPoly.resize(fNp);
    for (size_t i=0; i<a_n; ++i) {
      fPoly[i].X() = a_x[i];
      fPoly[i].Y() = a_y[i];
    }
    build_coeff(); // Build the spline coefficients.
  }

quintic() [3/3]

tools::spline::quintic::quintic ( const quintic &  a_from )

inline

Definition at line 436 of file spline.

:base_spline(a_from),fPoly(a_from.fPoly) {}

Member Function Documentation

build_coeff()

void tools::spline::quintic::build_coeff ( )

inlineprotected

algorithm 600, collected algorithms from acm. algorithm appeared in acm-trans. math. software, vol.9, no. 2, jun., 1983, p. 258-259.

quintic computes the coefficients of a quintic natural quintic spli s(x) with knots x(i) interpolating there to given function values: s(x(i)) = y(i) for i = 1,2, ..., n. in each interval (x(i),x(i+1)) the spline function s(xx) is a polynomial of fifth degree: s(xx) = ((((f(i)*p+e(i))*p+d(i))*p+c(i))*p+b(i))*p+y(i) (*) = ((((-f(i)*q+e(i+1))*q-d(i+1))*q+c(i+1))*q-b(i+1))*q+y(i+1) where p = xx - x(i) and q = x(i+1) - xx. (note the first subscript in the second expression.) the different polynomials are pieced together so that s(x) and its derivatives up to s"" are continuous.

input:

n number of data points, (at least three, i.e. n > 2) x(1:n) the strictly increasing or decreasing sequence of knots. the spacing must be such that the fifth power of x(i+1) - x(i) can be formed without overflow or underflow of exponents. y(1:n) the prescribed function values at the knots.

output:

b,c,d,e,f the computed spline coefficients as in (*). (1:n) specifically b(i) = s'(x(i)), c(i) = s"(x(i))/2, d(i) = s"'(x(i))/6, e(i) = s""(x(i))/24, f(i) = s""'(x(i))/120. f(n) is neither used nor altered. the five arrays b,c,d,e,f must always be distinct.

option:

it is possible to specify values for the first and second derivatives of the spline function at arbitrarily many knots. this is done by relaxing the requirement that the sequence of knots be strictly increasing or decreasing. specifically:

if x(j) = x(j+1) then s(x(j)) = y(j) and s'(x(j)) = y(j+1), if x(j) = x(j+1) = x(j+2) then in addition s"(x(j)) = y(j+2).

note that s""(x) is discontinuous at a double knot and, in addition, s"'(x) is discontinuous at a triple knot. the subroutine assigns y(i) to y(i+1) in these cases and also to y(i+2) at a triple knot. the representation (*) remains valid in each open interval (x(i),x(i+1)). at a double knot, x(j) = x(j+1), the output coefficients have the following values: y(j) = s(x(j)) = y(j+1) b(j) = s'(x(j)) = b(j+1) c(j) = s"(x(j))/2 = c(j+1) d(j) = s"'(x(j))/6 = d(j+1) e(j) = s""(x(j)-0)/24 e(j+1) = s""(x(j)+0)/24 f(j) = s""'(x(j)-0)/120 f(j+1) = s""'(x(j)+0)/120 at a triple knot, x(j) = x(j+1) = x(j+2), the output coefficients have the following values: y(j) = s(x(j)) = y(j+1) = y(j+2) b(j) = s'(x(j)) = b(j+1) = b(j+2) c(j) = s"(x(j))/2 = c(j+1) = c(j+2) d(j) = s"'((x(j)-0)/6 d(j+1) = 0 d(j+2) = s"'(x(j)+0)/6 e(j) = s""(x(j)-0)/24 e(j+1) = 0 e(j+2) = s""(x(j)+0)/24 f(j) = s""'(x(j)-0)/120 f(j+1) = 0 f(j+2) = s""'(x(j)+0)/120

Definition at line 479 of file spline.

                     {
 
    size_t i, m;
    double pqqr, p, q, r, s, t, u, v,
       b1, p2, p3, q2, q3, r2, pq, pr, qr;
 
    if (fNp <= 2) return;
 
    //     coefficients of a positive definite, pentadiagonal matrix,
    //     stored in D, E, F from 1 to n-3.
    m = fNp-2;
    q = fPoly[1].X()-fPoly[0].X();
    r = fPoly[2].X()-fPoly[1].X();
    q2 = q*q;
    r2 = r*r;
    qr = q+r;
    fPoly[0].D() = fPoly[0].E() = 0;
    if (q) fPoly[1].D() = q*6.*q2/(qr*qr);
    else fPoly[1].D() = 0;
 
    if (m > 1) {
       for (i = 1; i < m; ++i) {
          p = q;
          q = r;
          r = fPoly[i+2].X()-fPoly[i+1].X();
          p2 = q2;
          q2 = r2;
          r2 = r*r;
          pq = qr;
          qr = q+r;
          if (q) {
             q3 = q2*q;
             pr = p*r;
             pqqr = pq*qr;
             fPoly[i+1].D() = q3*6./(qr*qr);
             fPoly[i].D() += (q+q)*(pr*15.*pr+(p+r)*q
                                  *(pr* 20.+q2*7.)+q2*
                                  ((p2+r2)*8.+pr*21.+q2+q2))/(pqqr*pqqr);
             fPoly[i-1].D() += q3*6./(pq*pq);
             fPoly[i].E() = q2*(p*qr+pq*3.*(qr+r+r))/(pqqr*qr);
             fPoly[i-1].E() += q2*(r*pq+qr*3.*(pq+p+p))/(pqqr*pq);
             fPoly[i-1].F() = q3/pqqr;
          } else
             fPoly[i+1].D() = fPoly[i].E() = fPoly[i-1].F() = 0;
       }
    }
    if (r) fPoly[m-1].D() += r*6.*r2/(qr*qr);
 
    //     First and second order divided differences of the given function
    //     values, stored in b from 2 to n and in c from 3 to n
    //     respectively. care is taken of double and triple knots.
    for (i = 1; i < fNp; ++i) {
       if (fPoly[i].X() != fPoly[i-1].X()) {
          fPoly[i].B() =
             (fPoly[i].Y()-fPoly[i-1].Y())/(fPoly[i].X()-fPoly[i-1].X());
       } else {
          fPoly[i].B() = fPoly[i].Y();
          fPoly[i].Y() = fPoly[i-1].Y();
       }
    }
    for (i = 2; i < fNp; ++i) {
       if (fPoly[i].X() != fPoly[i-2].X()) {
          fPoly[i].C() =
             (fPoly[i].B()-fPoly[i-1].B())/(fPoly[i].X()-fPoly[i-2].X());
       } else {
          fPoly[i].C() = fPoly[i].B()*.5;
          fPoly[i].B() = fPoly[i-1].B();
       }
    }
 
    //     Solve the linear system with c(i+2) - c(i+1) as right-hand side.
    if (m > 1) {
       p=fPoly[0].C()=fPoly[m-1].E()=fPoly[0].F()
          =fPoly[m-2].F()=fPoly[m-1].F()=0;
       fPoly[1].C() = fPoly[3].C()-fPoly[2].C();
       fPoly[1].D() = 1./fPoly[1].D();
 
       if (m > 2) {
          for (i = 2; i < m; ++i) {
             q = fPoly[i-1].D()*fPoly[i-1].E();
             fPoly[i].D() = 1./(fPoly[i].D()-p*fPoly[i-2].F()-q*fPoly[i-1].E());
             fPoly[i].E() -= q*fPoly[i-1].F();
             fPoly[i].C() = fPoly[i+2].C()-fPoly[i+1].C()-p*fPoly[i-2].C()
                            -q*fPoly[i-1].C();
             p = fPoly[i-1].D()*fPoly[i-1].F();
          }
       }
    }
 
    fPoly[fNp-2].C() = fPoly[fNp-1].C() = 0;
    if (fNp > 3)
       for (i=fNp-3; i > 0; --i)
          fPoly[i].C() = (fPoly[i].C()-fPoly[i].E()*fPoly[i+1].C()
                         -fPoly[i].F()*fPoly[i+2].C())*fPoly[i].D();
 
    //     Integrate the third derivative of s(x)
    m = fNp-1;
    q = fPoly[1].X()-fPoly[0].X();
    r = fPoly[2].X()-fPoly[1].X();
    b1 = fPoly[1].B();
    q3 = q*q*q;
    qr = q+r;
    if (qr) {
       v = fPoly[1].C()/qr;
       t = v;
    } else
       v = t = 0;
    if (q) fPoly[0].F() = v/q;
    else fPoly[0].F() = 0;
    for (i = 1; i < m; ++i) {
       p = q;
       q = r;
       if (i != m-1) r = fPoly[i+2].X()-fPoly[i+1].X();
       else r = 0;
       p3 = q3;
       q3 = q*q*q;
       pq = qr;
       qr = q+r;
       s = t;
       if (qr) t = (fPoly[i+1].C()-fPoly[i].C())/qr;
       else t = 0;
       u = v;
       v = t-s;
       if (pq) {
          fPoly[i].F() = fPoly[i-1].F();
          if (q) fPoly[i].F() = v/q;
          fPoly[i].E() = s*5.;
          fPoly[i].D() = (fPoly[i].C()-q*s)*10;
          fPoly[i].C() =
          fPoly[i].D()*(p-q)+(fPoly[i+1].B()-fPoly[i].B()+(u-fPoly[i].E())*
                             p3-(v+fPoly[i].E())*q3)/pq;
          fPoly[i].B() = (p*(fPoly[i+1].B()-v*q3)+q*(fPoly[i].B()-u*p3))/pq-p
          *q*(fPoly[i].D()+fPoly[i].E()*(q-p));
       } else {
          fPoly[i].C() = fPoly[i-1].C();
          fPoly[i].D() = fPoly[i].E() = fPoly[i].F() = 0;
       }
    }
 
    //     End points x(1) and x(n)
    p = fPoly[1].X()-fPoly[0].X();
    s = fPoly[0].F()*p*p*p;
    fPoly[0].E() = fPoly[0].D() = 0;
    fPoly[0].C() = fPoly[1].C()-s*10;
    fPoly[0].B() = b1-(fPoly[0].C()+s)*p;
 
    q = fPoly[fNp-1].X()-fPoly[fNp-2].X();
    t = fPoly[fNp-2].F()*q*q*q;
    fPoly[fNp-1].E() = fPoly[fNp-1].D() = 0;
    fPoly[fNp-1].C() = fPoly[fNp-2].C()+t*10;
    fPoly[fNp-1].B() += (fPoly[fNp-1].C()-t)*q;
  }

eval()

double tools::spline::quintic::eval ( double  x ) const

inline

Definition at line 444 of file spline.

{if(!fNp) return 0;size_t klow=find_x(x);return fPoly[klow].eval(x);}

find_x()

size_t tools::spline::quintic::find_x ( double  x ) const

inlineprotected

Definition at line 447 of file spline.

                                {
    int klow=0;
 
    // If out of boundaries, extrapolate
    // It may be badly wrong
    if(x<=fXmin) klow=0;
    else if(x>=fXmax) klow=fNp-1;
    else {
      if(fKstep) { // Equidistant knots, use histogramming :
        klow = mn<int>(int((x-fXmin)/fDelta),fNp-1);
      } else {
        int khig=fNp-1;
        int khalf;
        // Non equidistant knots, binary search
        while((khig-klow)>1) {
          khalf = (klow+khig)/2;
          if(x>fPoly[khalf].X()) klow=khalf;
          else                   khig=khalf;
        }
      }
 
      // This could be removed, sanity check
      if( (x<fPoly[klow].X()) || (fPoly[klow+1].X()<x) ) {
        m_out << "tools::spline::quintic::find_x : Binary search failed"
              << " x(" << klow << ") = " << fPoly[klow].X() << " < x= " << x
              << " < x(" << klow+1<< ") = " << fPoly[klow+1].X() << "."
          << std::endl;
      }
    }
    return klow;
  }

operator=()

quintic& tools::spline::quintic::operator= ( const quintic &  a_from )

inline

Definition at line 437 of file spline.

                                            {
    if(this==&a_from) return *this;
    base_spline::operator=(a_from);
    fPoly = a_from.fPoly;
    return *this;
  }

Member Data Documentation

fPoly

std::vector<quintic_poly> tools::spline::quintic::fPoly

protected

Definition at line 697 of file spline.

---

The documentation for this class was generated from the following file:

• /Users/barrand/private/dev/softinex/g4tools/g4tools/tools/spline