How to integrate between 2 points of a given array.

22 次查看(过去 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
Please do not answer your own question, JUST to make a comment to add information. Worse, you then accepted that answer to your question. And finally, your "answer" did not add any information, since you had that same information in your question.
Answer by @Vahram Voskerchyan moved to a comment:
Problem with the indices should be
x2 = -0.4678
x1 = -0.4587
% find indices of these points
idx1 = find(field_5(:,1) >= x2, 1);
idx2 = find(field_5(:,1) >= x1, 1);
a1 = trapz(field_5(idx1:idx2,1),field_5(idx1:idx2,2)

请先登录,再进行评论。

采纳的回答

John D'Errico
John D'Errico 2021-11-27
编辑:John D'Errico 2021-11-27
idx1
idx1 =
2629
idx2
idx2 =
2622
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
idx1:idx2
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.
min(diff(field5(:,1)))
ans =
-3.7000
max(diff(field5(:,1)))
ans =
0.0100
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 =
2622
2623
2624
2625
2626
2627
2628
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.
numel(xs)
ans =
6026
numel(unique(xs))
ans =
2163
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 =
2163
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);
fnplt(spl,[x2,x1])
hold on
idx = find((xsc >= x2) & (xsc <= x1)); % plotting only points found in that interval
plot(xsc(idx),ysc(idx),'ro')
xline(x2,'g')
xline(x1,'g')
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 =
0.0146
Yet, if we use ONLY the data from those points inside x2 and x1, we would get:
trapz(xsc(idx),ysc(idx))
ans =
0.0117
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 =
0.0116
Find consolidator here:

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by