Hi,

I'm trying to do a cubic spline on the function y = (x-5)^2 with x = 0:9 and I'm having a problem linking the function file to the interp_1d.h and nr3.h files. Below is the code I've been using but it doesn't work when I pass it to g++. Any pointers for what I'm doing wrong would be really appreciated (I'm still quite new to C++)!

Thanks!

#include <iostream>
#include <iomanip>
#include <cmath>
#include "nr3.h"
#include "interp_1d.h"
using namespace std;

// Driver for routine interp_1d.h

int main() {
int xv [10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
int yv [10];
int ypn;
int yp1;
int i;
for (i = 1 ; i < 10 ; i++)
{
yv[i] = (xv[i] - 5)*(xv[i] - 5);
yp1 = -10;
ypn = 8;
}
yv=NR3::interp_1d(x, y, yp1, ypn);
}

...
I'm trying to do a cubic spline on the function y = (x-5)^2 with x = 0:9
Here's the drill:

Create NR3 VecDoubs with point values of x and y.


Use these (and values of yp1 and ypn if they are available) as arguments to construct a Spline_interp object.


Then to calculate the interpolated value of y, just feed a value of x to the interp() method of the Spline_interp object.


Like this:

//
// xspline_interp.cpp
// Example of Numerical Recipes version 3 spline interpolation
//
// davekw7x
//

#include "../code/nr3.h" // Use your path name to the NR3 distribution
#include "../code/interp_1d.h" // Use your path name to the NR3 distribution

// This is the function to be interpolated
double exact(double x)
{
return (x-5.0)*(x-5.0);
}

// Analytic derivative for calculating yp1 and ypn
double exactd(double x)
{
return 2*(x-5.0);
}

int main()
{
//
// Since Spline_interp needs vecDoubs we will
// use double precision arrays so that they
// can be "converted" by the vecDoub constructor.
//
double xv[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
const int num = sizeof(xv)/sizeof(xv[0]);
double yv[num];

for (int i = 0; i < num; i++)
{
yv[i] = exact(xv[i]);
}

// Maybe use these; maybe not
double yp1 = exactd(xv[0]);
double ypn = exactd(xv[num-1]);

// Now you have two arrays, xv and yv.
// The interpolation routine needs vecDoubs
VecDoub x(num, xv);
VecDoub y(num, yv);

// To use without yp1 and ypn, here's the call
//Spline_interp si = Spline_interp(x, y);

// Use yp1, ypn
Spline_interp si = Spline_interp(x, y, yp1, ypn);


// Will print n values of y = exact(x) and interpolated
// value of y for xmin <= x <= xmax
double xmin = 0.0;
double xmax = 9.0;
int n = 20;
double deltax = (xmax-xmin)/(n-1.0);


cout << " x exact interpolated" << endl;
cout << "--------------------------------" << endl;
cout << fixed;
for (int i = 0; i < n; i++) {
double xx = xmin + i*deltax;
cout << xx << " " << setw(10) << exact(xx) << " " << setw(10) << si.interp(xx) << endl;
}

cout << endl << "Taa-Daa!" << endl << endl;

return 0;

}


Output



x exact interpolated
--------------------------------
0.000000 25.000000 25.000000
0.473684 20.487535 20.487535
0.947368 16.423823 16.423823
1.421053 12.808864 12.808864
1.894737 9.642659 9.642659
2.368421 6.925208 6.925208
2.842105 4.656510 4.656510
3.315789 2.836565 2.836565
3.789474 1.465374 1.465374
4.263158 0.542936 0.542936
4.736842 0.069252 0.069252
5.210526 0.044321 0.044321
5.684211 0.468144 0.468144
6.157895 1.340720 1.340720
6.631579 2.662050 2.662050
7.105263 4.432133 4.432133
7.578947 6.650970 6.650970
8.052632 9.318560 9.318560
8.526316 12.434903 12.434903
9.000000 16.000000 16.000000

Taa-Daa!



Tested with GNU g++ version 4.1.2 on my Centos 5.8 Linux workstation. 23 October, 2012

Command line compilation with:

g++ -Wall -W xspline_interp.cpp -o xspline_interp



Regards,

Dave