What is the best way to go about integrating a discrete and evenly spaced dataset?

I am calculating a gas diffusional profile in a solid sphere for which the first and second derivatives of the curve are both negative. This causes a systematic bias such that linear interpolation/trapz causes volume calculation errors to rise asymptotically under certain circumstances. I have been told that Romberg integration is a good antidote to the problem, however, it appears that a requirement of the Romberg method is to have a smooth, continuous function. Unfortunately, this curve can not be described by such a function.

I am wondering if someone more skilled than myself could help me come up with an interpolation scheme (piecewise?) that would give me the continuous function that could be fed into a Romberg script.

I have also tried using Simpson's rule integration, but it seems my implementation is causing significant errors as well. I am probably doing something stupid wrong, but I have been staring at this stuff too long to really even tell. The approach here is to calculate and add the volume contributed by the concentration in each sphere between two spatial nodes, where v is the height of the concentration curve at radial position r, j is the index for spatial nodes and n is the number of spatial nodes:

    vol = simpsons([v(1) v(1)],0,r(1),[]) * (4/3 * pi * r(1)^3);
    for j = 2:n
        vol = vol + (4/3 * pi * r(j)^3 * simpsons(v(1:j),0,r(j),[]) - vol);
    end

I am probably making a stupid mistake here, and I'm not too familiar with interpolation schemes in MATLAB, so any advice or help on either front is greatly appreciated.

I have also attached an example of one of the curves I am trying to integrate over. v is on the y-axis and r (fractional) is on the x-axis.

Answer 1

[John D'Errico]

You have a relation with a derivative singularity. Newton-cotes panels, thus trapezoidal rules, simpson rules, etc. are all implicitly running a polynomial through the points, then integrating the polynomial. I recall showing that you can derive higher order Newton-Cotes rules from the lower order rules by such an extrapolation step. So Romberg won't really help much. It just generates an implicit higher order panel.

Rule 1: polynomials abhor singularities.

One option that I would ordinarily consider is to use an interpolating spline variant, such as pchip. Then integrate the pship interpolant. You can use fnint to do the integration. Do NOT use spline here. Use pchip. Classic splines are terrible with singularities of any sort, since they will generate nasty oscillations (ringing behavior.) Again, splines are based on polynomials, and what does rule 1 say?

You have not told us what order the singularity is. Since you know much about this relationship, you surely can know that behavior. Given that information, it would seem possible to devise a scheme that would work reasonably well. It would also be useful if you posted some data.

[Adam Goldsmith]

I'm not sure I'm following... most of what you are saying. What do you mean by singularity in this context?

The data itself includes tens of thousands of spatial nodes over millions of time steps, each with a slightly different curve (the picture above is a single timeslice). I can certainly post some examples, but I'm not sure what an efficient way to present such a large dataset here would be.

Also, I am reading the help files for pchip and fnint, but I quickly get lost as I don't have a background in this stuff. Would you mind explaining a bit of what these functions are doing, what their outputs are, and some examples of how they could be implemented?

edit: if by singularity you are referring to the steepness of the curve as the fractional radial position approaches 1, it only appears that way because it is simply a very, very steep slope. It becomes less and less steep with increasing time.

[Cedric]

Here are a few examples: https://blogs.mathworks.com/cleve/2012/07/16/splines-and-pchips/

[Adam Goldsmith]

Thanks, that's helpful but I still don't understand how to work with the struct output of pchip itself. How do I squeeze an actual piecewise function out of what it gives me for input into some integration scheme?

[Adam Goldsmith]

Still very confused as to how to extract a symbolic function from pchip to feed into some integration scheme... any help appreciated

[Cedric]

I cannot answer for John, and I didn't really read your question, but after scheming his answer I thought about something like:

    % - Build some data set and plot.
    x = [0, 0.5, 5, 10] ;       % Un-evenly distrib. => oscil. in spline interp.
    y = 1 + sin( x ) .* exp( -x/10 ) ;
    plot( x, y, 'kx', 'MarkerSize', 10 ) ;
    grid( 'on' ) ;
    hold( 'on' ) ;

    % - Build x for interp.
    xq = 0 : 0.2 : 10 ;

    % - Build spline interp. and plot.
    pp_spline = spline( x, y ) ;
    yq_spline = ppval( pp_spline, xq ) ;
    plot( xq, yq_spline, 'r-' ) ;

    % - Build pchip interp. and plot.
    pp_pchip  = pchip( x, y ) ;
    yq_pchip  = ppval( pp_pchip, xq ) ;
    plot( xq, yq_pchip, 'b-' ) ;

    legend( {'Data set', 'Spline interp.', 'Pchip interp.'} ) ;

    % - Integrate both.
    int_y_spline = integral( @(xx) ppval( pp_spline, xx), x(1), x(end) ) ;
    int_y_pchip  = integral( @(xx) ppval( pp_pchip, xx), x(1), x(end) ) ;
    fprintf( 'Int: spline = %.2f vs pchip = %2f.\n', int_y_spline, int_y_pchip ) ;

which outputs:

 Int: spline = 5.79 vs pchip = 7.421623.

where the piece of information that you were missing seems to be this PPVAL function.

[Cedric]

Actually, this was an old example and I updated it slightly to illustrate:

    % - Build some data set and plot.
    x = [0, 4, 4.9, 5.1, 10] ;
    y = 2 * (x < 5) - 1 ;
    subplot( 2, 2, 1 ) ;  grid( 'on' ) ;  hold( 'on' ) ;
    plot( x, y, 'kx', 'MarkerSize', 10 ) ;
    plot( x, y, 'k--', 'LineWidth', 1 ) ; 
    title( 'Data set' ) ;

    % - Build x for interp.
    xq = 0 : 0.1 : 10 ;

    % - Build spline interp. and plot.
    pp_spline = spline( x, y ) ;
    yq_spline = ppval( pp_spline, xq ) ;
    subplot( 2, 2, 2 ) ;  grid( 'on' ) ;  hold( 'on' ) ;
    plot( x, y, 'kx', 'MarkerSize', 10 ) ;
    plot( xq, yq_spline, 'r-' ) ;
    title( 'Spline interp.' ) ;

    % - Build pchip interp. and plot.
    pp_pchip  = pchip( x, y ) ;
    yq_pchip  = ppval( pp_pchip, xq ) ;
    subplot( 2, 2, 3 ) ;  grid( 'on' ) ;  hold( 'on' ) ;
    plot( x, y, 'kx', 'MarkerSize', 10 ) ;
    plot( xq, yq_pchip, 'b-' ) ;
    title( 'Pchip interp.' ) ;

    % - All together.
    subplot( 2, 2, 4 ) ;  grid( 'on' ) ;  hold( 'on' ) ;
    plot( x, y, 'kx', 'MarkerSize', 10 ) ;
    plot( xq, yq_spline, 'r-' ) ;
    plot( xq, yq_pchip, 'b-' ) ;
    legend( {'Data set', 'Spline interp.', 'Pchip interp.'}, ...
        'Location', 'southwest' ) ;

    % - Integrate both.
    int_y_spline = integral( @(xx) ppval( pp_spline, xx), x(1), x(end) ) ;
    int_y_pchip  = integral( @(xx) ppval( pp_pchip, xx), x(1), x(end) ) ;
    fprintf( 'Int: spline = %.2f vs pchip = %2f.\n', int_y_spline, int_y_pchip ) ;

The data is a sample of square signal in [-1, 1] whose integral over [0, 10] should be zero. Here is what we get:

 Int: spline = -215.79 vs pchip = 0.000000.

[Cedric]

or

 >> diff( fnval( fnint( pp_spline ), [0, 10] ))
 ans =
     -215.7866

 >> diff( fnval( fnint( pp_pchip ), [0, 10] ))
 ans =
     0