How to calculate are under the curve from 2 vector sets.
3 ビュー (過去 30 日間)
John D'Errico 2023 年 6 月 2 日
編集済み: John D'Errico 2023 年 6 月 2 日
It seems as if you want to compute the area between the curve and the x axis, while ignoring the sign of y. But a simple integral will cancel the positive and the negative parts. So essentially you want to compute the integral of the absolute value of your function.
There are easy ways to do this, and less easy. Do you want it easy, or good?
For example, since we don't have your data, consider this simple function:
F = @sin;
Now, I know the integral of @sin in one of those lobes is 2, so the integral of the absolute value will be 4 over the total interval. But if I just compute the integral on that total interval, any tool will give me zero, because the positive and negative lobes exactly cancel each other out.
But suppose all I have is a vector of data points?
x = linspace(0,2*pi,10);
y = F(x);
What you want is to compute the integral of abs(F). For example, if we have the function, integral can do it.
Indeed, integral was correct. (Oh, my god, so was I!)
However, if we try this:
then any integration rule that uses the absolute value on data points will probably poorly estimate the integral. If the spacing is fine enough, then it may not hurt a lot. But the integral there will be a bit too high because the region right around x==pi is screwed up. The funny thing is, trapz will tend to under-estimate the rest of integral, since it does that on functions that are negatively curved as is this one. So we have two competing sources of error.
intest = trapz(x,abs(y))
Not terrible, but not terribly accurate either. Approximately 1% error.
100*(intest - 4)/4
A better solution is to find the point where F crosses the x axis. Then interpolate your data points to find where that happens, and break the problem down into two integrals. And then finally, use a better integration rule than trapezoidal rule. (Even Simpson's wule would be good here.) The result would then be actually a pretty good estimate of the overall area.
Again, do you want a good solution, or do you want an easy solution? The above scheme would take some effort. Not really hugely much effort in my eyes. But some. Effort is a relative thing.
Ok, there is one trick we might try, that actually is pretty easy. First, interpolate your curve.
spl = spline(x,y);
xhat = linspace(min(x),max(x),1000);
yhat = ppval(spl,xhat);
intest2 = trapz(xhat,abs(yhat))
100*(intest2 - 4)/4
So this time, only 0.04 % error. The interpolation has reduced the error by a factor of roughly 25, and the use of spline here was trivially easy to do. The result might not be quite as good as what I might have gotten by use of a far more diligent code, but this was really quite easy.