Why does cumtrapz give lower magnitude for very thin Gaussian?

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?

Answer 1

[the cyclist]

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.

Answer 2

[John D'Errico]

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])

ans = 

1

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)

Ffun = function_handle with value:
    @(X,s)(sqrt(2.0).*1.0./sqrt(pi).*exp(X.^2.*1.0./s.^2.*(-1.0./2.0)))./(s.*2.0)

Ffun1 = @(x) Ffun(x,1)

Ffun1 = function_handle with value:
    @(x)Ffun(x,1)

Ffun01 = @(x) Ffun(x,0.01)

Ffun01 = function_handle with value:
    @(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)

ans = 
     0

trapz(x2,y2)

ans = 
          79.7884560802866

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.