How to integrate between 2 points of a given array.

4 次查看(过去 30 天)
I want to integrate between two points on the array, but when I do it I get 0.
x2 = -0.4678
x1 = -0.4587
% find indices of these points
idx1 = find(field5(:,1) >= x1, 1);
idx2 = find(field5(:,1) >= x2, 1);
a1 = trapz(field_5(idx1:idx2,1),field_5(idx1:idx2,2))
what I want to do is to integrate between calculated points of the graph.
  2 个评论



John D'Errico
John D'Errico 2021-11-27
编辑:John D'Errico 2021-11-27
idx1 =
idx2 =
So then what did you do? In order to understand (and then fix) the problem, you need to recognize what you are doing, and then to understand how MATLAB interptrets your actions.
a1 = trapz(field_5(idx1:idx2,1),field_5(idx1:idx2,2)
What is the vector
ans =
1×0 empty double row vector
That is because idx1 is GREATER than idx2. So the colon operator recognizes this must be an empty vector.
So you used trapz to integrate an empty array. What should it return? Logically, trapz might decide the result must be zero.
Should MATLAB know that perhaps you meant something different than the code you wrote? Should it be able to read your mind, to then decide that it should do differently, effectively rewriting your code? I might bet that if it made the wrong choice here when MATLAB decided to rewrite your code with no warning, your thoughts might validly include an expletive.
So how did idx1 and idx2 arise? To understand that, we need first to consider that field5(:,1) is NOT a sorted vector.
ans =
ans =
So merely finding the FIRST element of field5(:,1) that satisfies some property is not a good idea anyway.
Since I do not know what your real intention is here, I cannot give you a complete answer to your problem. Do you expect to integrate from x2 to x1, thus in a positive direction with increasing limits of integration? Are you hoping to integrate the elements from that variable from x2 to x1?
Finally, I might add that simply using trapz to perform that integration is asking for a terribly poor approximation, because if x1 and x2 are not in fact EXACTLY data points in that array, then there is part of the interval you are not using for integration. Far better MIGHT be to integrate a spline interpolant over that interval. I might do this:
idx = find((field5(:,1) >= x2) & (field5(:,1) <= x1))
idx =
Recognizing that this interval does not actually include x1 and x2, I might choose idx different. Thus, find the indicies from the largest point that does not exceed x2, then the smallest point that is not less than x1. And, since field5 is not a sorted array on the first column, better is first to sort it.
[xs,tags] = sort(field5(:,1));
ys = field5(tags,2);
This still suffers, because there are replicate elements in xs. While in the vicinity of x2 and x3, there are no reps, but even so, there are a LOT of replicates.
ans =
ans =
Best might be to average the reps down. Thus, replace each replicate x with a single element, with an average value for y at those points. I've posted the consolidator function on the file exchange for this purpose.
[xsc,ysc] = consolidator(xs,ys,@mean);
>> numel(xsc)
ans =
Finally, now we can build a spline interpolant of that curve. Since the curve appears to be a highly spiky thing, best is to use pchip.
spl = pchip(xsc,ysc);
hold on
idx = find((xsc >= x2) & (xsc <= x1)); % plotting only points found in that interval
grid on
You should see the problem now. If we integrate the function over the interval given only the indexes you found, using trapz, your integral will miss a fairly significant portion of the interval. So the true integral over that interval seems to be:
integral(@(x) ppval(spl,x),x2,x1)
ans =
Yet, if we use ONLY the data from those points inside x2 and x1, we would get:
ans =
So trapz was in error, underestimating the area by roughly 26%. Integral essentially agrees with trapz, if we restrict the interval. So the difference lies not in the use of a trapezoidal integration versus a spline, but in the use of the wrong interval.
integral(@(x) ppval(spl,x),xsc(idx(1)),xsc(idx(end)))
ans =
Find consolidator here:

更多回答(0 个)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by