Hi, I need help, I'm sure somebody has already solved similar problem. Here is the problem: I have file with data points (x, y) which represent function y = f(x) and i must calculate by numeriacal integration area under the curve.. So I know how to read data from file, i know how to use functions 'qtrap' or 'qsimp', i know how the interpolation works, but I don't know how to connect all these methods together..I have idea, that I will do interpolation thru data points and pass f(x) values thru 'func' to function 'qtrap', am I right? What are your recommendations? If you need more info, I can give it to you. Thak you..

...I have file with data points (x, y) which represent function y = f(x) and i must calculate by numeriacal integration area under the curve.. ...i know how to use functions 'qtrap' or 'qsimp', i know how the interpolation works...
I assume that you have read Chapters 3 and 4 of Numerical Recipes.

The functions qtrap and qsimp are composite methods that operate on points from a given function a step at a time, using smaller and smaller subintervals until the results from one trial to the next are sufficiently close.

If the underlying process that gave rise to the data points is anything close to a polynomial of degree 1 or 2, then, rather than using qtrap or qsimp, the Trapezoidal rule or Simpson's rule can be applied on the point set directly. (If the points are not equally spaced, Simpson's rule becomes a little more interesting, but still...) Depending on the nature of the data set, somewhat higher order Newton-Cotes formulas might be effectively employed (Simpson's 3/8 rule for example.)

Now, if you want to approximate the underlying process by some function induced from the data set, it seems to me that success might depend on what you know about the underlying process.

For example...

If you have, say, 20 data points, then because of numerical sensitivity, the 19th degree polynomial that passes exactly through the data points is (probably) a horrible approximation of the actual function, even if you knew that the underlying process is a 19th degree polynomial (and that, in itself, is very unlikely).

Sometimes people use a low-degree polynomial that satisfies some minimum-least-squares criterion for closeness to the data points. They plot out the approximating function and 'eyeball' it with respect to the data points and make a value judgment. (The raw number of the sum of the squares of the distances between data points and approximated points kind of leaves me cold, but it might be useful in some relative sense.)

On the other hand...

Polynomials might or might not be good candidates for approximating functions underlying physical processes. (Richard W. Hamming devoted more than a few pages in chapters of several books on the subject. In fact, there's a whole discipline based on successful use of trigonometric approximations.)

Anyhow...

Sometimes we assume that a physical system that gave rise to the data points is, somehow, smooth. We also might (have to) assume that the time between samples is short enough to capture any fast-changing data.

The concept of having a "smooth" function implies (to me) some kind of spline rather than a single polynomial that passes through (or near) all of the data points. (Unless, that is, we know that the underlying process is a low-degree polynomial.)

Bottom line: The key (to me) is not how to use qtrap or qsimp, but to understand enough about the underlying process to get reasonable results. (And that implies that you have some way of deciding what is "reasonable.") See Footnote.

Regards,

Dave


Footnote:

I didn't answer your question, did I?

Well, the steps that I see would be...

Decide on some approximating function (polynomial or not) for your data set. This has to be the most important step, right?

Use some method to get an estimation of the parameters of the approximating function. You might find one of the methods from the "Modeling of Data" chapter of Numerical Recipes to be appropriate.

Create a program function that uses these parameters to evaluate the function at a given data point.

Feed this function along with beginning and ending data points to qsimp or qtrap.

If this doesn't make sense (or if it, somehow, it doesn't work), then post again with some more specific information.


Maybe someone can help above and beyond the vague generalities of my post (but there are no guarantees). In particular, tell us about what version of Numerical Recipes you are using. How are you compiling? Show us the program that you are using. Attach a text file with data points. Stuff like that.

Thak you for your answer.
I'm using polynomial interpolation thru 4 points or rational interp. thru 4 points and it looks like a good choice. I've written some code that gives me some acceptable result, but I don't know if it's right. Here is the code:

/* Assigment
* 1. [Numerical integration, interpolation] Function y = f(x) is given by table in file "func.dat" (left column 'x' and right column 'y').
* Calculate by numerical integration value of integral from 0 to 1 f(x) dx. Actually discrete values given by the table are calculated
* by substituting into analytically specified function, which definite integral thru interval <0, 1> is easy to calculate and which is with
* precision at 9 points valid to 0.502451028. Compare this exact value with results obtained by using various methods of interpolation and
* numerical integration.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NRANSI
#include "nrutil.h"
#include "nr.h"
#define N 86 //number of points
#define step 128
#define file "func.txt"

float *xa, *ya, y, dy;
long unsigned int v;
int k, n;

float func(float x) {

float h = 1.0 / step;
int i;

for (i = 0; i <= step; i++) { // Interpolation

x = i * h;
hunt(xa, N, x, &v);
k = IMIN(IMAX(v - (n - 1) / 2, 1), N + 1 - n);
ratint(xa + k - 1, ya + k - 1, n, x, &y, &dy);

}

return y/2.0;

}

int main() {

int i, j;
float a = 0.0, b = 1.0, s;
FILE *input;

printf("Insert number of points for interpolation: ");
scanf("%d",&n);

xa = vector(1, N);
ya = vector(1, N);

if ((input = fopen(file,"r")) == NULL)
nrerror("Data file not found\n");


for (j = 1; j <= 86; j++) { // Reading from file

fscanf(input, "%f %f", &xa[j], &ya[j]);

}

fclose(input);

s = qtrap(func, a, b);

printf("Value of integral using qtrap: %.9f\n", s);
printf("Given value of integral is 0.502451028.\n");

free_vector(xa, 1, N);
free_vector(ya, 1, N);

return 0;
}
#undef NRANSI
I'm using 2nd edition of Numerical recipes. Compiler is gcc on latest Xubuntu, I'm supporter of classical Makefile compiling of code. Data points are in the attached file 'func.txt'.

...some code that gives me some acceptable result, but I don't know...

Hmmm... I can't exactly see the point of your code.

Try something like the following:

Change your func to print the value it returns each time it is called:

float func(float x)
{
float xx = x;
float h = 1.0 / step;
int i;

for (i = 0; i <= step; i++) {
x = i * h;
hunt(xa, N, x, &v);
k = IMIN(IMAX(v - (n - 1) / 2, 1), N + 1 - n);
ratint(xa + k - 1, ya + k - 1, n, x, &y, &dy);
}

printf("func(%f) returning %f\n", xx, y/2.0);
return y/2.0;
}



Here's the output I got for a run:

Insert number of points for interpolation: 4
func(0.000000) returning 0.502073
func(1.000000) returning 0.502073
func(0.500000) returning 0.502073
func(0.250000) returning 0.502073
func(0.750000) returning 0.502073
func(0.125000) returning 0.502073
func(0.375000) returning 0.502073
func(0.625000) returning 0.502073
func(0.875000) returning 0.502073
func(0.062500) returning 0.502073
func(0.187500) returning 0.502073
func(0.312500) returning 0.502073
func(0.437500) returning 0.502073
func(0.562500) returning 0.502073
func(0.687500) returning 0.502073
func(0.812500) returning 0.502073
func(0.937500) returning 0.502073
func(0.031250) returning 0.502073
func(0.093750) returning 0.502073
func(0.156250) returning 0.502073
func(0.218750) returning 0.502073
func(0.281250) returning 0.502073
func(0.343750) returning 0.502073
func(0.406250) returning 0.502073
func(0.468750) returning 0.502073
func(0.531250) returning 0.502073
func(0.593750) returning 0.502073
func(0.656250) returning 0.502073
func(0.718750) returning 0.502073
func(0.781250) returning 0.502073
func(0.843750) returning 0.502073
func(0.906250) returning 0.502073
func(0.968750) returning 0.502073
Value of integral using qtrap: 0.502073169
Given value of integral is 0.502451028.


Now, whether the final result is something you are happy with or not, it is obvious to me that qtrap is not contributing anything here. The function func is supposed to be an approximating function that returns a the (approximate) value of y for a given value of x. Yours doesn't do it for me.

I made a cubic spline approximating function and used at as follows:

#include <stdio.h>
#include <stdlib.h>

#include "nrutil.h"
#include "nr.h"

/*
These need to be global, since they are used
in main() and also in func(). There is no way
to pass them as parameters to func(), since
qtrap and qsimp require func to have a single
argument.
*/
const int NPOINTS = 86;
float *xa, *ya, *y2a;

float func(float);

int main()
{

int i;
float a = 0.0, b = 1.0, integralValue;
FILE *infile;
const char *inname = "func.dat";

printf("Given value of integral = %f\n", 0.502451028);

xa = vector(1, NPOINTS);
ya = vector(1, NPOINTS);
y2a = vector(1, NPOINTS);

if ((infile = fopen(inname,"r")) == NULL) {
printf("Can't open file %s for reading.\n", inname);
return EXIT_FAILURE;
}

for (i = 1; i <= NPOINTS; i++) {
if (fscanf(infile, "%f %f", &xa[i], &ya[i]) != 2) {
printf("Problem reading point %d\n", i);
return EXIT_FAILURE;
}
}
fclose(infile);

/*
Fill the array y2a so that splint can return interpolated
values.
*/
spline(xa, ya, NPOINTS, 1.0e31, 1.0e31, y2a);

integralValue = qtrap(func, a, b);
printf("Value of integral using qtrap = %f\n", integralValue);

integralValue = qsimp(func, a, b);
printf("Value of integral using qsimp = %f\n", integralValue);

free_vector(xa, 1, NPOINTS);
free_vector(ya, 1, NPOINTS);
free_vector(y2a, 1, NPOINTS);

return 0;
}

float func(float x)
{
float y;
/* Use previously calculated array y2a to get an
interpolated value of y for a given x
*/
splint(xa, ya, y2a, NPOINTS, x, &y);

/* Uncomment the following for debugging */
/* printf("func(%f) returning %f\n", x, y); */
return y;
}

Output on my Linux system (compiled with GNU gcc version 4.1.2):

Given value of integral = 0.502451
Value of integral using qtrap = 0.502445
Value of integral using qsimp = 0.502444


Regards,

Dave

Footnote:
Since I have no insight as to the actual underlying process, I don't know whether a polynomial or a rational polynomial function is appropriate. If I did, I might try to find "best approximation" of the parameters with one of the modeling techniques from Chapter 15 of Numerical Recipes.

Lacking actual mathematical justification to try polynomial stuff, I just used a cubic spline approximation to get a continuous function through data points that is "smooth" everywhere on the interval of integration. This may or may not be analytically superior, but at least it shows a rather straightforward way to get results from qtrap or qsimp. (Note that equivalent functions from the NR3 distribution use double precision arithmetic, and both qtrap and qsimp give values of 0.502444 to six significant digits.)

Thank you very much Dave. Now I can see, that my code didn't make sense. Your code and advices showed me the way how to solve similar problems. But I still need to read Chapter 15 from NR, to I knew something about Data Modeling.
with kind regards,
Marcel

Here is my solution using rational interpolation:

#include <stdio.h>
#include <stdlib.h>

#include "nrutil.h"
#include "nr.h"

const int NP = 86, n = 6;
int k;
long unsigned int v;
float *xa, *ya, dy;

float func(float);

int main() {

int i;
float a = 0.0, b = 1.0, integralValue;
FILE *input;
const char *fname = "func.dat";

printf("Given value of integral = 0.502451028\n");

xa = vector(1, NP);
ya = vector(1, NP);

if ((input = fopen(fname,"r")) == NULL) {
printf("Can't open file %s for reading.\n", fname);
return EXIT_FAILURE;
}

for (i = 1; i <= NP; i++) {
if (fscanf(input, "%f %f", &xa[i], &ya[i]) != 2) {
printf("Problem reading point %d\n", i);
return EXIT_FAILURE;
}
}

fclose(input);


integralValue = qtrap(func, a, b);
printf("Value of integral using qtrap = %.9f\n", integralValue);

integralValue = qsimp(func, a, b);
printf("Value of integral using qsimp = %.9f\n", integralValue);

integralValue = qromb(func, a, b);
printf("Value of integral using qromb = %.9f\n", integralValue);

free_vector(xa, 1, NP);
free_vector(ya, 1, NP);

return 0;
}

float func(float x) {

float y;

hunt(xa, NP, x, &v);
k = IMIN(IMAX(v - (n - 1) / 2, 1), NP + 1 - n);
ratint(xa + k, ya + k, n, x, &y, &dy);
//printf("x = %f \t y = %f\n", x, y);

return y;

}