本文共 8221 字，大约阅读时间需要 27 分钟。

/* * spline.h * * simple cubic spline interpolation library without external * dependencies * * --------------------------------------------------------------------- * Copyright (C) 2011, 2014 Tino Kluge (ttk448 at gmail.com) * *  This program is free software; you can redistribute it and/or *  modify it under the terms of the GNU General Public License *  as published by the Free Software Foundation; either version 2 *  of the License, or (at your option) any later version. * *  This program is distributed in the hope that it will be useful, *  but WITHOUT ANY WARRANTY; without even the implied warranty of *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the *  GNU General Public License for more details. * *  You should have received a copy of the GNU General Public License *  along with this program.  If not, see 
   . * --------------------------------------------------------------------- * */#ifndef TK_SPLINE_H#define TK_SPLINE_H#include 
   
    #include 
    
     #include 
     
      #include 
      
       // unnamed namespace only because the implementation is in this// header file and we don't want to export symbols to the obj filesnamespace{namespace tk{// band matrix solverclass band_matrix{private:    std::vector< std::vector
       
         > m_upper;  // upper band    std::vector< std::vector
        
          > m_lower; // lower bandpublic: band_matrix() {}; // constructor band_matrix(int dim, int n_u, int n_l); // constructor ~band_matrix() {}; // destructor void resize(int dim, int n_u, int n_l); // init with dim,n_u,n_l int dim() const; // matrix dimension int num_upper() const { return m_upper.size()-1; } int num_lower() const { return m_lower.size()-1; } // access operator double & operator () (int i, int j); // write double operator () (int i, int j) const; // read // we can store an additional diogonal (in m_lower) double& saved_diag(int i); double saved_diag(int i) const; void lu_decompose(); std::vector
         
           r_solve(const std::vector
          
           & b) const; std::vector
           
             l_solve(const std::vector
            
             & b) const; std::vector
             
               lu_solve(const std::vector
              
               & b, bool is_lu_decomposed=false);};// spline interpolationclass spline{public: enum bd_type { first_deriv = 1, second_deriv = 2 };private: std::vector
               
                 m_x,m_y; // x,y coordinates of points // interpolation parameters // f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i std::vector
                
                  m_a,m_b,m_c; // spline coefficients double m_b0, m_c0; // for left extrapol bd_type m_left, m_right; double m_left_value, m_right_value; bool m_force_linear_extrapolation;public: // set default boundary condition to be zero curvature at both ends spline(): m_left(second_deriv), m_right(second_deriv), m_left_value(0.0), m_right_value(0.0), m_force_linear_extrapolation(false) { ; } // optional, but if called it has to come be before set_points() void set_boundary(bd_type left, double left_value, bd_type right, double right_value, bool force_linear_extrapolation=false); void set_points(const std::vector
                 
                  & x, const std::vector
                  
                   & y, bool cubic_spline=true); double operator() (double x) const; double deriv(int order, double x) const;};// ---------------------------------------------------------------------// implementation part, which could be separated into a cpp file// ---------------------------------------------------------------------// band_matrix implementation// -------------------------band_matrix::band_matrix(int dim, int n_u, int n_l){ resize(dim, n_u, n_l);}void band_matrix::resize(int dim, int n_u, int n_l){ assert(dim>0); assert(n_u>=0); assert(n_l>=0); m_upper.resize(n_u+1); m_lower.resize(n_l+1); for(size_t i=0; i
                   
                    0) { return m_upper[0].size(); } else { return 0; }}// defines the new operator (), so that we can access the elements// by A(i,j), index going from i=0,...,dim()-1double & band_matrix::operator () (int i, int j){ int k=j-i; // what band is the entry assert( (i>=0) && (i
                    
                     =0) && (j
                     
                       diogonal, k<0 lower left part, k>0 upper right part if(k>=0) return m_upper[k][i]; else return m_lower[-k][i];}double band_matrix::operator () (int i, int j) const{ int k=j-i; // what band is the entry assert( (i>=0) && (i
                      
                       =0) && (j
                       
                         diogonal, k<0 lower left part, k>0 upper right part if(k>=0) return m_upper[k][i]; else return m_lower[-k][i];}// second diag (used in LU decomposition), saved in m_lowerdouble band_matrix::saved_diag(int i) const{ assert( (i>=0) && (i
                        
                         =0) && (i
                         
                          dim(); i++) { assert(this->operator()(i,i)!=0.0); this->saved_diag(i)=1.0/this->operator()(i,i); j_min=std::max(0,i-this->num_lower()); j_max=std::min(this->dim()-1,i+this->num_upper()); for(int j=j_min; j<=j_max; j++) { this->operator()(i,j) *= this->saved_diag(i); } this->operator()(i,i)=1.0; // prevents rounding errors } // Gauss LR-Decomposition for(int k=0; k
                          
                           dim(); k++) { i_max=std::min(this->dim()-1,k+this->num_lower()); // num_lower not a mistake! for(int i=k+1; i<=i_max; i++) { assert(this->operator()(k,k)!=0.0); x=-this->operator()(i,k)/this->operator()(k,k); this->operator()(i,k)=-x; // assembly part of L j_max=std::min(this->dim()-1,k+this->num_upper()); for(int j=k+1; j<=j_max; j++) { // assembly part of R this->operator()(i,j)=this->operator()(i,j)+x*this->operator()(k,j); } } }}// solves Ly=bstd::vector
                           
                             band_matrix::l_solve(const std::vector
                            
                             & b) const{ assert( this->dim()==(int)b.size() ); std::vector
                             
                               x(this->dim()); int j_start; double sum; for(int i=0; i
                              
                               dim(); i++) { sum=0; j_start=std::max(0,i-this->num_lower()); for(int j=j_start; j
                               
                                operator()(i,j)*x[j]; x[i]=(b[i]*this->saved_diag(i)) - sum; } return x;}// solves Rx=ystd::vector
                                
                                  band_matrix::r_solve(const std::vector
                                 
                                  & b) const{ assert( this->dim()==(int)b.size() ); std::vector
                                  
                                    x(this->dim()); int j_stop; double sum; for(int i=this->dim()-1; i>=0; i--) { sum=0; j_stop=std::min(this->dim()-1,i+this->num_upper()); for(int j=i+1; j<=j_stop; j++) sum += this->operator()(i,j)*x[j]; x[i]=( b[i] - sum ) / this->operator()(i,i); } return x;}std::vector
                                   
                                     band_matrix::lu_solve(const std::vector
                                    
                                     & b, bool is_lu_decomposed){ assert( this->dim()==(int)b.size() ); std::vector
                                     
                                       x,y; if(is_lu_decomposed==false) { this->lu_decompose(); } y=this->l_solve(b); x=this->r_solve(y); return x;}// spline implementation// -----------------------void spline::set_boundary(spline::bd_type left, double left_value, spline::bd_type right, double right_value, bool force_linear_extrapolation){ assert(m_x.size()==0); // set_points() must not have happened yet m_left=left; m_right=right; m_left_value=left_value; m_right_value=right_value; m_force_linear_extrapolation=force_linear_extrapolation;}void spline::set_points(const std::vector
                                      
                                       & x, const std::vector
                                       
                                        & y, bool cubic_spline){ assert(x.size()==y.size()); assert(x.size()>2); m_x=x; m_y=y; int n=x.size(); // TODO: maybe sort x and y, rather than returning an error for(int i=0; i
                                        
                                          rhs(n); for(int i=1; i
                                         
                                          ::const_iterator it; it=std::lower_bound(m_x.begin(),m_x.end(),x); int idx=std::max( int(it-m_x.begin())-1, 0); double h=x-m_x[idx]; double interpol; if(x
                                          
                                           m_x[n-1]) { // extrapolation to the right interpol=(m_b[n-1]*h + m_c[n-1])*h + m_y[n-1]; } else { // interpolation interpol=((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx]; } return interpol;}double spline::deriv(int order, double x) const{ assert(order>0); size_t n=m_x.size(); // find the closest point m_x[idx] < x, idx=0 even if x
                                           
                                            ::const_iterator it; it=std::lower_bound(m_x.begin(),m_x.end(),x); int idx=std::max( int(it-m_x.begin())-1, 0); double h=x-m_x[idx]; double interpol; if(x
                                            
                                             m_x[n-1]) { // extrapolation to the right switch(order) { case 1: interpol=2.0*m_b[n-1]*h + m_c[n-1]; break; case 2: interpol=2.0*m_b[n-1]; break; default: interpol=0.0; break; } } else { // interpolation switch(order) { case 1: interpol=(3.0*m_a[idx]*h + 2.0*m_b[idx])*h + m_c[idx]; break; case 2: interpol=6.0*m_a[idx]*h + 2.0*m_b[idx]; break; case 3: interpol=6.0*m_a[idx]; break; default: interpol=0.0; break; } } return interpol;}} // namespace tk} // namespace#endif /* TK_SPLINE_H */
                                            
                                           
                                          
                                         
                                        
                                       
                                      
                                     
                                    
                                   
                                  
                                 
                                
                               
                              
                             
                            
                           
                          
                         
                        
                       
                      
                     
                    
                   
                  
                 
                
               
              
             
            
           
          
         
        
       
      
     
    
   

这个是开源项目（https://github.com/ttk592/spline）的代码。这里声明一下。

下面将使用：

 std::vector<double> X,Y;

X->pushback(数据);

 Y->pushback(数据);

。。。。。。。。。。。。。。

tk::spline s;  //定义一个函数。

s.set_points(X,Y); //将已知的点集合放入函数。

cout<<x<<s(x); //x就是你累加的坐标x，自己使用循环增加，而s(x)就是该x轴值对应的y轴的坐标了。

转载地址：https://islet.blog.csdn.net/article/details/78498209 如侵犯您的版权，请留言回复原文章的地址，我们会给您删除此文章，给您带来不便请您谅解！