splines.cpp

#include <vector>
#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

class Spline {
    public:
        vector<double> x, f, c;
        Spline(const string & filename);
        double calc(double x1, double & p, double & p1, double & p2) const;
};

int main() {
    Spline s("04_in.txt");
    cout << "x\tf(x)\ts(x)\teps" << endl;
    for (int i = 0, n = s.x.size()-1; i < n; ++i) {
        double p, p1, p2;
        s.calc(s.x[i], p, p1, p2);
        double eps = fabs(s.f[i]-p)/fabs(s.f[i]);
        cout << s.x[i] << "\t" << s.f[i] << "\t" << p << "\t" << eps << endl;
    }
    double p, p1, p2;
    s.calc(1.5, p, p1, p2);
    cout << "s(1.5) = " << p << endl;
    cout << "f(1.5*pi/6) = " << cos(1.5*M_PI/6) << endl;
    return 0;
}
// f(i, x) = a(i)+b(i)*(x-x[i-1])+c(i)*(x-x[i-1])^2+d(i)*(x-x[i-1])^3
Spline::Spline(const string & filename) {
    ifstream fin(filename);
    int n;
    fin >> n;
    x.assign(n+1, 0);
    f.assign(n+1, 0);
    c.assign(n+1, 0);
    for (int i = 0; i < n; ++i)
        fin >> x[i] >> f[i];

    vector<double> k(n+1, 0);
    k[1] = c[1] = 0;
    for (int i = 2; i <= n; ++i) {
        int j = i - 1, m = j - 1;
        double a, b, r;
        a = x[i] - x[j];
        b = x[j] - x[m];
        r = 2*(a+b)-b*c[j];
        c[i] = a/r;
        k[i] = (3.0*((f[i]-f[j])/a-(f[j]-f[m])/b)-b*k[j])/r;
    }

    c[n] = k[n];
    for (int i = n-1; i >= 2; --i)
        c[i] = k[i]-c[i]*c[i+1];
}

double Spline::calc(double x1, double & p, double & p1, double & p2) const {
    int i = 1, n = x.size();
    while (x1 > x[i] && i != n) ++i;
    int j = i-1;
    double a, b, q, r, d;


    printf("\t--> i=%d:\t", i);
    printf("a=%lf\t", f[i-1]);
    printf("b=%lf\t", (f[i]-f[i-1])/(x[i]-x[i-1])-(c[i+1]+2*c[i])*(x[i]-x[i-1])/3.0);
    printf("c=%lf\t", c[i]);
    printf("d=%lf\t\n", (c[i+1]-c[i])/(3.0*(x[i]-x[i-1])));
    a = f[j];
    b = x[j];
    q = x[i] - b;
    r = x1-b;
    p = c[i];
    d = c[i+1];

    b = (f[i]-a)/q-(d+2*p)*q/3.0;
    d = (d-p)/q*r;


    p1 = b + r*(2*p+d);
    p2 = 2*(p+d);
    p = a+r*(b+r*(p+d/3.0));
}