Is dynamic interpolation of input of function for 'quad' integration possible?

My task is to integrate a function given limited data, e.g.

max_energy = 1000;  
quad(@(x) a.*x.*exp(-b.*x),0,max_energy)

where a, b are vectors, functions of energy. My problem is that it seems quad uses Simpson's rule iteratively until some tolerance is met, so the array size of a, b would need to change dynamically to match. Is it possible to interpolate between the data points I have to match quad? Must I do it in advance or use a different integration method?

Answer 1

[Andrew Newell]

Given that your data are irregularly spaced, you may have to use interpolation after all. If your points are (xdata,ydata), you could integrate as below:

pp = interp1(xdata,ydata,'pchip','pp');
f = @(x) ppval(pp,x);
max_energy = 1000;  
quad(f,0,max_energy)

Just don't assume the integral is as accurate as the default tolerance.

[Daniel]

I'm not certain it's appropriate to use Hermite polynomials for interpolation, rather than linearly (the table lookup method, right?) ... I'm trying the latter, I think, simply using interp1(xdata,ydata,integration_variable) within quad and am getting NaN (not a number?) as the answer. Is there some overflow issue occurring since I'm having it interpolate a value for every point of integration? The code:

muinterp = @(xdata,ydata,int_var) interp1(xdata,ydata,int_var);

quad(@(x) muinterp(energy,murho,x),0,1)

The result:

Warning: Infinite or Not-a-Number function value encountered.

> In quad at 109

ans = NaN

[Daniel]

The method you've provided works; my 'default' interp1 use does not. Does it use linear interpolation if you don't specify a method? Do you know what problem the method I tried has that the 'pchip' method doesn't have?

[Andrew Newell]

The default method is linear. I don't know why you're getting the NaN, but if you use interp in the form yi = interp1(x,Y,xi), you have to keep recalculating the interpolating polynomials.

Answer 2

[Andrew Newell]

Any increase in accuracy based on interpolating a and b would probably be illusory. If your points are regularly spaced, you could use trapz:

y = a.*x.*exp(-b.*x);
yInt = trapz(x,y);

[Bård Skaflestad]

That's a good point. If the 'a' and 'b' functions are known only empirically, then |trapz| is the way to go. Don't forget to scale the result with the sample point distance if the integrand is not sampled at unit intervals.

[Andrew Newell]

Actually, scaling is only necessary if you use yInt = trapz(y). Consider this example:

x = 1:.1:10;

y = x.^2;

yInt = cumtrapz(x,y);

plot(x,yInt,x,x.^3/3)

The documentation is confusing, though - I had to read it a couple of times before I realized that scaling is not necessary.

[Daniel]

My data is not evenly spaced unless I divide it into segments and work the problem piecemeal; cf. http://physics.nist.gov/PhysRefData/XrayMassCoef/ComTab/water.html . That seems difficult to do without losing information or discarding certain data points. Would interpolating to get evenly spaced data be useful? But it seems doing so would make trapz no more accurate than quad.

[Bård Skaflestad]

@Andrew

Thanks a lot -- I'd missed that part of the documentation. I obviously need to pay more attention...

[Andrew Newell]

I send MATLAB some feedback about the documentation.

Answer 3

[Bård Skaflestad]

All of MATLAB's quadrature methods require an integrand that can be evaluated at vector inputs and return an equally sized vector result.

Do you mean to say that your a and b in some way depend on x? If so, you may have to implement a traditional function (i.e., .m) file that evaluates both a and b along with the resulting integrand.

For instance,

function y = integrand(x)
  a = some_function(x);
  b = some_other_function(x);
  y = a .* x .* exp(-b .* x);
end

Does this help at all?

[Bård Skaflestad]

And then, of course, I forgot the |quad| call:

quad(@integrand, 0, max_energy)

[Daniel]

It does help, thank you, but I don't think I have an expression for a(x), b(x) -- yes, their values depend on the variable of integration -- that instead their values are experimentally determined, hence the need for interpolation. Is curve fitting the way to go, then, based on the existing data for a, b? But the data is seemingly logarithmic -- will fitting to polynomials work for something whose behavior changes significantly depending on the range of x? (and I have never done curve fitting before, so I would need to learn; would you recommend resources?)