What is the best way to go about integrating a discrete and evenly spaced dataset?
7 次查看(过去 30 天)
显示 更早的评论
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/167462/image.png)
0 个评论
回答(1 个)
John D'Errico
2017-9-23
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.
7 个评论
Cedric
2017-9-24
编辑:Cedric
2017-9-24
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/185463/image.png)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spline Postprocessing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!