// 文件：cbspline-test.cpp
// C++ Cubic Spline Interpolation using the Boost library.
// Interpolation is an estimation of a value within two known values in a sequence of values.
// From a selected test function y(x) = (5.0/(x))*sin(5*x)+(x-6), 21 numbers of (x,y) equal-spaced sequential points were generated.
// Using the known 21 (x,y) points, 1001 numbers of (x,y) equal-spaced cubic spline interpolated points were calculated.
// Since the test function y(x) is known exactly, the exact y-value for any x-value can be calculated.
// For each point, the interpolation error is calculated as the difference between the exact (x,y) value and the interpolated (x,y) value.
// This program writes outputs as tab delimited results to both screen and named text file
// COMPILATION: $ g++ -lgmp -lm -std=c++11 -o cbspline-test.xx cbspline-test.cpp
// EXECUTION: $ ./cbspline-test.xx
// GNUPLOT 1: gnuplot> load "cbspline-versus-exact-function-plots.gpl"
// GNUPLOT 2: gnuplot> load "cbspline-interpolated-point-errors.gpl"
#include <iostream>
#include <fstream>
#include <boost/math/interpolators/cubic_b_spline.hpp>
using namespace std;
// ========================================================
int main(int argc, char* argv[]) {
// ========================================================
// Vector size 21 = 20 space intervals, size 1001 = 1000 space intervals
std::vector<double> x_knot(21), x_exact(1001), x_cbspline(1001);
std::vector<double> y_knot(21), y_exact(1001), y_cbspline(1001);
std::vector<double> y_diff(1001);
double x; double x0 = 0.0; double t0 = 0.0;
double xstep1 = 0.5; double xstep2 = 0.01;
ofstream outfile; // Output data file tab delimited values
outfile.open("cbspline-errors-1000-points.dat"); // y_diff = (y_exact - y_cbspline)
// Handling zero-point infinity (1/x) when x = 0
x0 = 1e-18;
x_knot[0] = x0;
y_knot[0] = (5.0/(x0))*sin(5*x0)+(x0-6); // Selected test function
for (int i = 1; i <= 20; i++) { // Note: Loop i begins with 1 not 0, 20 points
x = (i*xstep1);
x_knot[i] = x;
y_knot[i] = (5.0/(x))*sin(5*x)+(x-6);
}
// Execution of the cubic spline interpolation from Boost library
// NOTE: Using xstep1 = 0.5 based on knot points stepsize (for 21 known points)
boost::math::cubic_b_spline<double> spline(y_knot.begin(), y_knot.end(), t0, xstep1);
// Again handling zero-point infinity (1/x) when x = 0
x_cbspline[0] = x_knot[0];
y_cbspline[0] = y_knot[0];
for (int i = 1; i <= 1000; ++i) { // 1000 points.
x = (i*xstep2);
x_cbspline[i] = x;
// NOTE: y = spline(x) is valid for any value of x in xrange [0:10]
// meaning, any x within range [x_knot.begin() and x_knot.end()]
y_cbspline[i] = spline(x);
}
// Write headers for screen display and output file
std::cout << "# x_exact[i] \t y_exact[i] \t y_cbspline[i] \t (y_diff[i])" << std::endl;
outfile << "# x_exact[i] \t y_exact[i] \t y_cbspline[i] \t (y_diff[i])" << std::endl;
// Again handling zero-point infinity (1/x) when x = 0
x_exact[0] = x_knot[0];
y_exact[0] = y_knot[0];
y_diff[0] = (y_exact[0] - y_cbspline[0]);
std::cout << x_exact[0] << "\t" << y_exact[0] << "\t" << y_cbspline[0] << "\t" << y_diff[0] << std::endl; // Write to screen
outfile << x_exact[0] << "\t" << y_exact[0] << "\t" << y_cbspline[0] << "\t" << y_diff[0] << std::endl; // Write to text file
for (int i = 1; i <= 1000; ++i) { // 1000 points
x = (i*xstep2);
x_exact[i] = x;
y_exact[i] = (5.0/(x))*sin(5*x)+(x-6);
y_diff[i] = (y_exact[i] - y_cbspline[i]);
std::cout << x_exact[i] << "\t" << y_exact[i] << "\t" << y_cbspline[i] << "\t" << y_diff[i] << std::endl; // Write to screen
outfile << x_exact[i] << "\t" << y_exact[i] << "\t" << y_cbspline[i] << "\t" << y_diff[i] << std::endl; // Write to text file
}
outfile.close();
return(0);
}
// ========================================================
/*
# GNUPLOT script 1
# File: cbspline-versus-exact-function-plots.gpl
set term qt persist size 700,500
set xrange [-1:10.0]
set yrange [-12.0:20.0]
# set autoscale # set xtics automatically
set xtics 0.5
set ytics 2.0
# set xtic auto # set xtics automatically
# set ytic auto # set ytics automatically
set grid x
set grid y
set xlabel "x"
set ylabel "y(x)"
set title "Function points y(x) = (5.0/(x))*sin(5*x)+(x-6)"
set yzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
set xzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
show xzeroaxis
show yzeroaxis
show grid
plot 'cbspline-errors-1000-points.dat' using 1:2 with lines linetype 3 linewidth 1.0 linecolor 'blue' title 'exact function curve', 'cbspline-errors-1000-points.dat' using 1:3 with lines linecolor 'red' linetype 3 linewidth 1.0 title 'cubic spline interpolated curve'
*/
// ========================================================
/*
# GNUPLOT script 2
# File: cbspline-interpolated-point-errors.gpl
set term qt persist size 700,500
set xrange [-1:10.0]
set yrange [-2.50:2.5]
# set autoscale # set xtics automatically
set xtics 0.5
set ytics 0.5
# set xtic auto # set xtics automatically
# set ytic auto # set ytics automatically
set grid x
set grid y
set xlabel "x"
set ylabel "y"
set title "Function points y = (5.0/(x))*sin(5*x)+(x-6)"
set yzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
set xzeroaxis linetype 3 linewidth 1.5 linecolor 'black'
show xzeroaxis
show yzeroaxis
show grid
plot 'cbspline-errors-1000-points.dat' using 1:4 with lines linetype 3 linewidth 1.0 linecolor 'red' title 'cubic spline interpolated errors'
*/
// ========================================================
