Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Are there any good libraries for solving cubic splines in C++? [closed]

Tags:

I'm looking for a good C++ library to give me functions to solve for large cubic splines (on the order of 1000 points) anyone know one?

like image 605
Faken Avatar asked Jul 30 '09 05:07

Faken


People also ask

Is cubic spline interpolation the best?

Runge's phenomenon tells us that such an approximation often has large oscillations near the ends of the interpolating interval. On the other hand, cubic spline interpolation is often considered a better approximation method because it is not prone to such oscillations.

Which Matlab command can be used for to construct natural cubic spline for the given data set?

s = spline( x , y , xq ) returns a vector of interpolated values s corresponding to the query points in xq . The values of s are determined by cubic spline interpolation of x and y .


1 Answers

Write your own. Here is spline() function I wrote based on excellent wiki algorithm:

#include<iostream>
#include<vector>
#include<algorithm>
#include<cmath>
using namespace std;

using vec = vector<double>;

struct SplineSet{
    double a;
    double b;
    double c;
    double d;
    double x;
};

vector<SplineSet> spline(vec &x, vec &y)
{
    int n = x.size()-1;
    vec a;
    a.insert(a.begin(), y.begin(), y.end());
    vec b(n);
    vec d(n);
    vec h;

    for(int i = 0; i < n; ++i)
        h.push_back(x[i+1]-x[i]);

    vec alpha;
    alpha.push_back(0);
    for(int i = 1; i < n; ++i)
        alpha.push_back( 3*(a[i+1]-a[i])/h[i] - 3*(a[i]-a[i-1])/h[i-1]  );

    vec c(n+1);
    vec l(n+1);
    vec mu(n+1);
    vec z(n+1);
    l[0] = 1;
    mu[0] = 0;
    z[0] = 0;

    for(int i = 1; i < n; ++i)
    {
        l[i] = 2 *(x[i+1]-x[i-1])-h[i-1]*mu[i-1];
        mu[i] = h[i]/l[i];
        z[i] = (alpha[i]-h[i-1]*z[i-1])/l[i];
    }

    l[n] = 1;
    z[n] = 0;
    c[n] = 0;

    for(int j = n-1; j >= 0; --j)
    {
        c[j] = z [j] - mu[j] * c[j+1];
        b[j] = (a[j+1]-a[j])/h[j]-h[j]*(c[j+1]+2*c[j])/3;
        d[j] = (c[j+1]-c[j])/3/h[j];
    }

    vector<SplineSet> output_set(n);
    for(int i = 0; i < n; ++i)
    {
        output_set[i].a = a[i];
        output_set[i].b = b[i];
        output_set[i].c = c[i];
        output_set[i].d = d[i];
        output_set[i].x = x[i];
    }
    return output_set;
}

int main()
{
    vec x(11);
    vec y(11);
    for(int i = 0; i < x.size(); ++i)
    {
        x[i] = i;
        y[i] = sin(i);
    }

    vector<SplineSet> cs = spline(x, y);
    for(int i = 0; i < cs.size(); ++i)
        cout << cs[i].d << "\t" << cs[i].c << "\t" << cs[i].b << "\t" << cs[i].a << endl;
}
like image 129
cpp Avatar answered Oct 20 '22 06:10

cpp