interpolation.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//interpolation.hpp 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 3 of the License, or
//(at your option) any later version.
//
//interpolation.hpp 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 interpolation.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
// Copyright(C) 2015 - 2021 by Bertini2 Development Team
//
// See <http://www.gnu.org/licenses/> for a copy of the license, 
// as well as COPYING.  Bertini2 is provided with permitted 
// additional terms in the b2/licenses/ directory.
 
// individual authors of this file include:
// silviana amethyst, university of wisconsin eau claire
// Tim Hodges, Colorado State University

 
#pragma once
 
#include "bertini2/common/config.hpp"
 
namespace bertini{
    namespace endgame{
          
template<typename CT>       
    Vec<CT> HermiteInterpolateAndSolve(CT const& target_time, const unsigned int num_sample_points, const TimeCont<CT> & times, const SampCont<CT> & samples, const SampCont<CT> & derivatives, ContStart shift_from = ContStart::Back)
{
    assert((times.size() >= num_sample_points) && "must have sufficient number of sample times");
    assert((samples.size() >= num_sample_points) && "must have sufficient number of sample points");
    assert((derivatives.size() >= num_sample_points) && "must have sufficient number of derivatives");
 
    
    unsigned num_t, num_s, num_d;
    if (shift_from == ContStart::Back)
    {
        num_t = times.size()-1;
        num_s = samples.size()-1;
        num_d = derivatives.size()-1;
    }
    else
    {
        num_t = num_s = num_d = num_sample_points-1;
    }
 
    Mat< Vec<CT> > space_differences(2*num_sample_points,2*num_sample_points);
    Vec<CT> time_differences(2*num_sample_points);
 
 
    for(unsigned int ii=0; ii<num_sample_points; ++ii)
    { 
        space_differences(2*ii,0)   = samples[    num_s-ii];        /*  F[2*i][0]    = samples[i];    */
        space_differences(2*ii+1,0) = samples[    num_s-ii];    /*  F[2*i+1][0]  = samples[i];    */
        space_differences(2*ii+1,1) = derivatives[num_d-ii];        /*  F[2*i+1][1]  = derivatives[i]; */
        time_differences(2*ii)      = times[      num_t-ii];        /*  z[2*i]       = times[i];       */
        time_differences(2*ii+1)    = times[      num_t-ii];        /*  z[2*i+1]     = times[i];       */
    }
 
    //Add first round of finite differences to fill out rest of matrix. 
    for(unsigned int ii=1; ii< num_sample_points; ++ii)
    {
        space_differences(2*ii,1) = (space_differences(2*ii,0) - space_differences(2*ii-1,0)) / (time_differences(2*ii) - time_differences(2*ii-1));
    
    }
 
    //Filling out finite difference matrix to get the diagonal for hermite interpolation polyonomial.
    for(unsigned int ii=2; ii < 2*num_sample_points; ++ii)
    {
        for(unsigned int jj=2; jj <=ii; ++jj)
        {
            space_differences(ii,jj) = 
                (space_differences(ii,jj-1) - space_differences(ii-1,jj-1)) 
                    / 
                (time_differences(ii) - time_differences(ii-jj));                       
        }
    }
 
    //Start of Result from Hermite polynomial, this is using the diagonal of the 
    //finite difference matrix.
    Vec<CT> Result = space_differences(2*num_sample_points - 1,2*num_sample_points - 1); 
 
 
    //This builds the hermite polynomial from the highest term down. 
    //As we multiply the previous result we will construct the highest term down to the last term.
    for (unsigned ii=num_sample_points-1; ii >= 1; --ii)
    {
        Result = ((Result*(target_time - time_differences(ii)) + space_differences(2*ii, 2*ii)) * (target_time - time_differences(ii-1)) + space_differences(2*ii-1, 2*ii-1)).eval();  
    }
    
    // Last term in hermite polynomial.
    return (Result * (target_time - time_differences(0)) + space_differences(0,0)).eval(); 
} //re: HermiteInterpolateAndSolve
 
}}  // re: namespaces