src/mathlib_src/GCubic_Spline_Interpolation.cpp

説明を見る。

#include "GCubic_Spline_Interpolation.h"
#include <iostream>

using namespace GMathLib;

GCubic_Spline_Interpolation::GCubic_Spline_Interpolation()
: GObject(), p_cfData(0), p_at(0), p_u(0), p_v(0), p_h(0)
{
    Class_Name("GCubic_Spline_Interpolation");
}

GCubic_Spline_Interpolation::~GCubic_Spline_Interpolation()
{
    if(p_cfData != 0){
        delete[] p_cfData;
        delete p_at;
        delete p_h;
        delete p_u;
        delete p_v;
    }
}

void GCubic_Spline_Interpolation::Prepair(int gridp_num)
{
    if(p_cfData != 0){
        delete[] p_cfData;
        delete p_at;
        delete p_h;
        delete p_u;
        delete p_v;
    }
    
    //gridp_num：格子点の数(N+1個)
    p_cfData = new GCubicFncData[gridp_num-1];
    p_at = new GMatrix(gridp_num, gridp_num); //行列A(N+1*N+1)
    p_u = new GVector(gridp_num, GVector::COLUMN_VECTOR); //列ベクトルu_i(N+1)
    p_v = new GVector(gridp_num, GVector::COLUMN_VECTOR); //列ベクトルv_i(N+1)
    p_h = new GVector(gridp_num-1, GVector::COLUMN_VECTOR); //h_i(N)
}

int GCubic_Spline_Interpolation::DoInterpolation(GMatrix& dc_mx)
{
    
    int d_size;

    if(p_cfData == 0){
        std::cout << "Failed to Interpolation. Call Prepair() before calling DoInterpolation()!!" << std::endl;

        return 1;
    }

    if(dc_mx.Column() != 2 || dc_mx.Row() != p_at->Row()){
        std::cout << "Failed to Interpolation. Invalid matrix size" << std::endl;
        return 1;
    }

    //d_sizeに列ベクトルの要素数(N+1)にセットする
    d_size = p_at->Column(); 

    //差分幅 h を計算し, 格子点間の幅情報を保持する
    for(int i=0; i < d_size-1; i++){
        p_h->Set(i, dc_mx(i+1, 0) - dc_mx(i, 0));       
    }

    // N+1*N+1の行列Aの設定

    //まず行列要素の端を設定
    p_at->Fill(0);
    p_at->Set(0, 0, 1.0);
    p_at->Set(d_size-1, d_size-1, 1.0);
    //残りの行列要素（対角要素）の設定
    for(int i=1; i < d_size-1; i++){
        p_at->Set(i, i-1, (*p_h)(i-1));
        p_at->Set(i, i, 2.0 * ( (*p_h)(i-1) + (*p_h)(i)) );
        p_at->Set(i, i+1, (*p_h)(i));
    }
    

    //N+1の列ベクトルCの設定
    //まず両端を設定
    p_v->Set(0, 0);
    p_v->Set(d_size-1, 0);
    //残りの行列要素を設定
    double tmp = 0.0;
    for(int i=1; i < d_size-1; i++){
        tmp =  3.0 * 
               (
                 ( dc_mx(i+1, 1) - dc_mx(i, 1)) / (*p_h)(i) - ( dc_mx(i, 1) - dc_mx(i-1, 1)) / (*p_h)(i-1)
               );
        p_v->Set(i, tmp);
    }   
    
    //Matrix(A) * Vector(u)=Vector(v)をGMathLib::GLinear_Homo_Eqクラスを使って，
    //列ベクトルvについて解く．
    //at->Print();
    ln_hm_solver.S_Gauss_Partial_Pivot(*p_at, *p_v, *p_u);


    //求めた列ベクトルuを使って，各区間の3次多項式の係数データを求め，
    //GCubicFncData配列に格納する
    for(int i=0; i < d_size-1; i++){
        (p_cfData + i)->a = ( (*p_u)(i+1) - (*p_u)(i) ) / ( 3.0 * (*p_h)(i));
        (p_cfData + i)->b = (*p_u)(i);       
        (p_cfData + i)->c = ( dc_mx(i+1, 1) - dc_mx(i,1) ) / (*p_h)(i) 
                            - 
                            (*p_h)(i) * ( 2.0 * (*p_u)(i) + (*p_u)(i+1) ) / 3.0;

        (p_cfData + i)->d = dc_mx(i, 1);
    }
      
    return 0;    
}

---