Why does cumtrapz give lower magnitude for very thin Gaussian?
조회 수: 6 (최근 30일)
이전 댓글 표시
I want to generate a cumulative integral of a Gaussian curve. Using a nice broad curve gives a magnitude I would expect:
clear
% t domain
t_min = 0;
t_max = 10;
numx = 1E3;
t_samples = t_max/(numx-1);
x = t_min:t_samples:t_max;
% Variables
scale = 1;
shift = t_max/4; %t_max/2 to centre
width = 1;
T_o = 0;
a = zeros(size(x));
for i = 1:length(x)
a(i) = scale * exp( (-(x(i) - shift) ^2)/width) + T_o;
end
plot(x,a)
y = cumtrapz(x,a);
plot(x,y)
But once I set up the Gaussian to have dimensions that we observe in experiment (1ps wide laser pulse with 46µJ of energy), cumtrapz returns a magnitude much lower than the height of the original Gaussian.
clear
% t domain
t_min = 0;
t_max = 1E-11;
numx = 1E3;
t_samples = t_max/(numx-1);
x = t_min:t_samples:t_max;
% Variables
scale = 46E-6;
shift = t_max/16; %t_max/2 to centre
width = 5E-26;
T_o = 0;
a = zeros(size(x));
for i = 1:length(x)
a(i) = scale * exp( (-(x(i) - shift) ^2)/width) + T_o;
end
plot(x,a)
y = cumtrapz(x,a);
plot(x,y)
I'm guessing this has something to do with trapezoidal integration not being precise, but I feel like I'm missing something. How do I go about making cumtrapz in the second case accurate?
댓글 수: 0
채택된 답변
the cyclist
2023년 4월 12일
편집: the cyclist
2023년 4월 12일
I see no surprising behavior here. Combining the plots into one, here are your results:
% t domain
t_min = 0;
t_max = 10;
numx = 1E3;
t_samples = t_max/(numx-1);
x = t_min:t_samples:t_max;
% Variables
scale = 1;
shift = t_max/4; %t_max/2 to centre
width = 1;
T_o = 0;
a = zeros(size(x));
for i = 1:length(x)
a(i) = scale * exp( (-(x(i) - shift) ^2)/width) + T_o;
end
figure
tiledlayout(2,2)
nexttile(1)
plot(x,a)
nexttile(3)
y = cumtrapz(x,a);
plot(x,y)
% t domain
t_min = 0;
t_max = 1E-11;
numx = 1E3;
t_samples = t_max/(numx-1);
x = t_min:t_samples:t_max;
% Variables
scale = 46E-6;
shift = t_max/16; %t_max/2 to centre
width = 5E-26;
T_o = 0;
a = zeros(size(x));
for i = 1:length(x)
a(i) = scale * exp( (-(x(i) - shift) ^2)/width) + T_o;
end
nexttile(2)
plot(x,a)
nexttile(4)
y = cumtrapz(x,a);
plot(x,y)
In both cases, cumtrapz is calculating the (cumulative) area under the curve. In your first case (left side of my figure), the magnitude of both your x values and y values are order 1, so the area being order 1 is expected. In your second case (right side of my figure), your x is of order 10^-12, and your y is of order 10^-5, so the area being of order 10^-17 is expected.
I don't know why you expected otherwise.
댓글 수: 0
추가 답변 (1개)
John D'Errico
2023년 4월 12일
편집: John D'Errico
2023년 4월 12일
As @the cyclist tells you. I would suggest you accept his answer as the correct one. I'm posting an answer here, as opposed to a comment, because comments often get lost and missed.
Remember that cumtrapz is NOT an exact integration tool. It computes a trapezoidal integration, on the set of points you gave it. If those points are too widely spaced to yield an accurate estimate of the integral, then it has a problem. In the example you have given that underestimates the area, a peak that is infinitely high, but infinitely narrow will be a difficult thing to integrate. Expect almost any adaptive integration tool to fail here. And that very possibly could include integral itself. (Int because it is a symbolic tool, can survive.)
For example
syms X real
syms s real positive
F = exp(-(X/s)^2/2)/s/sqrt(2*sym(pi));
int(F,X,[-inf,inf])
And we expect the integral to be exactly 1, as it is known to be, for any value of s. The problem is, compare what we would get for a trapezoidal integral when s is small. to that for larger values of s.
Ffun = matlabFunction(F)
Ffun1 = @(x) Ffun(x,1)
Ffun01 = @(x) Ffun(x,0.01)
fplot(Ffun1,[-5,5])
hold on
fplot(Ffun01,[-5,5])
hold off
The point is, both of those functions contain the same area under the curve. One is wide and flat, and easy to integrate accurately, the other is beginning to be a good approximation to a Dirac delta function.
Look more closely at the narrow one, for two sets of x values.
x1 = -5:2:5;
x2 = -6:2:6;
y1 = Ffun01(x1);
y2 = Ffun01(x2);
fplot(Ffun01,[-5,5])
hold on
plot(x1,y1,'ro-')
plot(x2,y2,'bx-')
hold off
I made the two sets of nodes quite coarse here, to make the point. The red curve skips over x==0, so it completely misses the entire mass of the function.
format long g
trapz(x1,y1)
trapz(x2,y2)
A trapezoidal rule (and this includes cumtrapz) does nothing more intelligently than compute the area under the connect-the-dots curve. If that connect-the-dots curve is a poor approximation to the function, then expect it to either wildly under-predict or over-predict the true integral.
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!