Difference between integral and cumtrapz.

I have a tempral data series, which I need to integrate. In order to understand how Matlab does that, I start with a parabola as the data series. Could you guys please explain me why I get completely different results using the integral and cumtrapz funtions? The cumtrapz is giving me a cubic function, but I have no idea why it is only positive and I also do not understand the limits. My understanding is that the answer should be simple: intfx = x^3/3...

x = [-10:0.01:10]';
fx = x.^2;
% intfx = integral(@(x) fx, x(1), x(end), 'ArrayValued', true);
intfx = cumtrapz(fx,x);
plot(x,fx,'b')
hold on
plot(x,intfx,'r')

Answer 1

[Star Strider]

They produce the same result if you let them.

Try something like this --

x = [-10:0.01:10]';
fx = @(x) x.^2;
for k = 1:numel(x)
    intfx_int(k) = integral(fx, x(1), x(k), 'ArrayValued', true);
end

intfx_ctz = cumtrapz(x,fx(x));

figure
tiledlayout(2,1)
nexttile
plot(x,fx(x),'b')
hold on
plot(x,intfx_ctz,'r')
hold off

nexttile
plot(x,fx(x),'b')
hold on
plot(x,intfx_int,'r')
hold off

●

Your original code will produce a single value for 'intfx' (that I use as 'intfx_int') since you give it two integration limits, so it will integrate between them to produce a single result that is the integral of the function between those limits:

intfx = integral(@(x)x.^2, x(1), x(end), 'ArrayValued', true)

intfx = 666.6667

Here, I gave it varying limits of integration in the loop.

.

[Aram Zeytunyan]

Thank you all for your help and recommendations. However, I would like to ask a follow-up question. What would be the right way to integrate an imported data? I guess in that case I cannot use the funtional form any more (fx = @(x) x.^2). The imported data is a measured spectrum, and I would like to integrate within a specific range of wavelengths. I have attached the data file.

[Star Strider]

You are integrating a vector, so use the cumtrapz function.

Perhaps something like this --

spectrum = readtable('spectrum.txt');
spectrum.Properties.VariableNames=["Frequency","Magnitude"]

spectrum = 3648×2 table
    Frequency    Magnitude
    _________    _________

     188.83        -35.5  
     188.89        -13.2  
     188.96        -37.1  
     189.03        -33.6  
     189.09          5.5  
     189.16         13.2  
     189.22        -16.7  
     189.29         16.2  
     189.36          2.6  
     189.42         33.3  
     189.49         -1.3  
     189.56         17.2  
     189.62         13.9  
     189.69         43.2  
     189.75         -7.7  
     189.82         31.2  

spectrumIntegral = cumtrapz(spectrum.Frequency, spectrum.Magnitude);

figure
yyaxis left
plot(spectrum.Frequency, spectrum.Magnitude)
xlabel('Frequency')
ylabel('Magnitude')
yyaxis right
plot(spectrum.Frequency, spectrumIntegral)
ylabel('Integrated Magnitude')
grid

●

spectrumIntegralSquared = cumtrapz(spectrum.Frequency, spectrum.Magnitude.^2);

figure
yyaxis left
plot(spectrum.Frequency, spectrum.Magnitude)
xlabel('Frequency')
ylabel('Magnitude')
yyaxis right
plot(spectrum.Frequency, spectrumIntegralSquared)
ylabel('Integrated Magnitude Squared')
grid

I assume the 'range of wavelengths' you want are the ones you have chosen top ut in the 'spectrum.txt' file. If you want to restrict them further, you can do that with 'logical indexing', that is relatively straightforward.

Make appropriate changes to get the result you want.

EDIT --

Corrected typographical errors.

.

[Aram Zeytunyan]

Thank you! I was hoping to use the integral funtion to be able to set the integration limits, but apparently that's not how it works. Maybe I should have started this thread with my real use case from the very beginning.

I need to calculate the integral from the attached image (Equation 3), where x is the measured temporal profile, f(nu) is the measured spectral profile, g0 and L are constants.

With no spectral dependence (Eq. 4), I am getting the result that I would expect. But when I try to calculate the integral, I am getting a curve similar to what you showed, which has no physical meaning.

This is the line in my code that properly models Eq. 4.

with_no_sp = (exp((g-alpha)*L)-1).*(g-g*alpha/g0)./(g-alpha);

But the integral returns a wrong result.

with_sp1 = abs(cumtrapz(freq, (g.*spectr.*(exp((g.*spectr-alpha)*L)-1))./((g.*spectr-alpha)*L)));

Strangely enough, when I put freq at the end, it gives me a result that is very close to the expectation...

with_sp2 = cumtrapz((g.*spectr.*(exp((g.*spectr-alpha)*L)-1))./((g.*spectr-alpha)*L),freq);

[Star Strider]

As always, my pleasure!

You can set the integration limits with 'logical indexing', although I am not certain what limits you want to impose here.

imshow(imread('integral.png'))

What you want to do seems to require that you first model the function by fitting it to the data (I'm not certain what the variables mean here, so I can't do that using the information I currently have).

My impression from the image is that you want to do some ort of Fourier transform (or perhaps inverse Fourier transform) on the data. To do an inverse Fourier transform would require that you supply a complex vector as the argument. How that vector is depicted (one-sided or two-sided Fourier transform) would then dictate how the rest of the code uses it and inverts it.

A Fourier transform of it (since there is only a real vector) would go something like this --

spectrum = readtable('spectrum.txt');
spectrum.Properties.VariableNames=["Frequency","Magnitude"]

spectrum = 3648×2 table
    Frequency    Magnitude
    _________    _________

     188.83        -35.5  
     188.89        -13.2  
     188.96        -37.1  
     189.03        -33.6  
     189.09          5.5  
     189.16         13.2  
     189.22        -16.7  
     189.29         16.2  
     189.36          2.6  
     189.42         33.3  
     189.49         -1.3  
     189.56         17.2  
     189.62         13.9  
     189.69         43.2  
     189.75         -7.7  
     189.82         31.2  

% spectrumIntegral = cumtrapz(spectrum.Frequency, spectrum.Magnitude);

[FTs1,Fv] = FFT1(spectrum.Frequency, spectrum.Magnitude);

figure
% yyaxis left
loglog(Fv, abs(FTs1))
xlabel('Time')
ylabel('Amplitude')
% yyaxis right
% plot(spectrum.Frequency, spectrumIntegral)
% ylabel('Integrated Magnitude')
grid
title('Fourier Transform of ''spectrum'' Data')

% 
% spectrumIntegralSquared = cumtrapz(spectrum.Frequency, spectrum.Magnitude.^2);
% 
% figure
% yyaxis left
% plot(spectrum.Frequency, spectrum.Magnitude)
% xlabel('Frequency')
% ylabel('Magnitude')
% yyaxis right
% plot(spectrum.Frequency, spectrumIntegralSquared)
% ylabel('Integrated Magnitude Squared')
% grid


function [FTs1,Fv] = FFT1(s,t)
% One-Sided Numerical Fourier Transform
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    Fv = Fv(:);
    FTs1 = FTs(Iv,:);
end

This is simply a guess on my part. I do not understand what you actually want to do. I also do not understand the variables in the image.

.

[Aram Zeytunyan]

The equations are describing light amplification in a gain medium, where L is the lenght of the medium, g0 is the gain, and alpha is the loss coefficient. I have done a lot of Fourier transforms between temporal and spectral profiles, but this is not the case. In this case, they are independent. The temporal profile of the amplified light needs to be calculated taking into account the spectral dependence of the gain.

My code is the following. Variables amp_temporal and amp_temporal_nosp should be Equations 3 and 4, respectively. As you can see, amp_temporal_nosp (which is equation 4) works fine, but the integral for amp_temporal is not returning any meaningful result (it should be similar to amp_temporal_nosp, just slightly asymmetric and with higher secondary peak). The physics is straightforward, it's just the integration that is not working.

clear
close all

timedata = readtable('temporal.csv');
temporal = timedata.Var2;
temporal = temporal/max(temporal);
tstart = -200.2; % ns, side
tstop = 199.8; % ns, side
t = [tstart:(tstop-tstart)/(length(temporal)-1):tstop]';

spectraldata = readtable('spectrum.txt');
lambda0 = spectraldata.Var1; % nm
spectr0 = spectraldata.Var2;
spectr0 = spectr0-min(spectr0);
spectr0 = spectr0/max(spectr0);
longcut = 197;
% longcut = 400;
shortcut = 180;
lambda000 = lambda0(lambda0>shortcut & lambda0<longcut);
spectr000 = spectr0(lambda0>shortcut & lambda0<longcut);
lambda0 = lambda000;
spectr0 = spectr000;
lambda = [lambda0(1):(lambda0(end)-lambda0(1))/(length(temporal)-1):lambda0(end)]';
freq = 2.998*10^8./lambda; % ns-1
% fx = fit(lambda0,spectr0,'sin8');
% spfitted = feval(fx,lambda0);
% fx = fit(freq0,spectr0,'sin8');
% spfitted = feval(fx,freq0);
% spfitted(spfitted < 0) = 0;
spectr = spline(lambda0,spectr0,lambda);
spectr = spectr/max(spectr);
% spectr = flip(spectr);

g0 = 0.14; % cm^-1;
L = 54; % cm
alpha = 7*0.008; % cm^-1
g = temporal*g0;
% amp_temporal = abs(cumtrapz((g.*spectr.*(exp((g.*spectr-alpha)*L)-1))./((g.*spectr-alpha)*L)));
amp_temporal = abs(cumtrapz(freq,(g.*spectr.*(exp((g.*spectr-alpha)*L)-1))./((g.*spectr-alpha)*L)));
amp_temporal_nosp = (exp((g-alpha)*L)-1).*(g-g*alpha/g0)./(g-alpha);

figure(1)
plot(t,temporal,'b','LineWidth',1);
xlabel('time (ns)')
ylabel('signal (mV)')
xlim([-50 100]) % ns
ylim([-0.01 1])
hold on
plot(t,amp_temporal/max(amp_temporal),'r','LineWidth',2);
plot(t,amp_temporal_nosp/max(amp_temporal_nosp),'k','LineWidth',2);

figure(2)
plot(lambda0,spectr0,'b','LineWidth',1);
hold on
plot(lambda,spectr,'r','LineWidth',1);
xlabel('wavelength (nm)')
ylabel('signal (counts)')
% xlim([189 413])
xlim([189 max(lambda)])
ylim([0 max(spectr)])

[Star Strider]

I admit to being a bit lost.

I note that in the .png, 'amp_temporal' (Equation 3) is the product of two separate integrations. However in your code (unless I'm missing something), it's calculated using only one integration.

I don't fully understand what you're doing with respect to how the variables map to the .png equations, (and I'm not familiar with the physics), so I won't attempt to code that myself here.

[Aram Zeytunyan]

The second integral is a constant, there is no x in it. If we believe that equation 3 is correct, there is no need to write a long line for the denominator in the equation, it's not going to change the shape. Thanks again for your help. I now know that cumtrapz is the right/only way to calculate the integral.

[Star Strider]

As always, my pleasure!

It doesn't look like a constant to me since it integrates over ν (and it multiplies as the inverse, so essentially the first integral is divided by it), however I don't understand either the physics or the variable definitions, so I'll defer to you.

Answer 2

[Walter Roberson]

integral() takes a function handle and does an adaptive quadrature numeric integration of the given function. The adaptive quadrature is exact for polynomials up to roughly degree 32767 if I read the source code correctly. For non-polynomials, the accuracy will depend on how close the approximating polynomial gets to the actual function.

cumtrapz() uses trapezoid numeric approximation, which is only accurate for polynomials up to (I think) degree 4 that are sampled often enough.

So, cumtrapz() will be less precise than integral() for all but very low degree polynomials.

The sample function you are integrating is a very low degree polynomial so the results should be the same to within numeric round-off.

Answer 3

[Torsten]

Your input arguments to "cumtrapz" are reversed: you have to use

intfx = cumtrapz(x,fx);

instead of

intfx = cumtrapz(fx,x);

And the answer with "cumtrapz" as well as with "integral" will be x^3/3 - (-10)^3/3 because both start with 0 at x = -10.

[Star Strider]

I reversed them o be correct inmy answer. I did not mention the change.

[John D'Errico]

This has always bothered me about trapz and cumtrapz. That is, if you call them with one argument, it is assumed to be y, and the stride in x is assumed to be 1. So we can do this:

x = 1:10;
y = x.^3 - x + 1;;
cumtrapz(y)

ans = 1×10

0 4 20 63 154 320 594 1015 1628 2484

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now, if I add additional information such as a non-unit stride between the x values, then logically it should be with a second argument, which if not supplied, has a default value. So I would expect this to be the case, if I want to tell cumtrapz both the vectors x and y.

cumtrapz(y,x)

ans = 1×10

0 9 54 180 450 945 1764 3024 4860 7425

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

cumtrapz(x,y)

ans = 1×10

0 4 20 63 154 320 594 1015 1628 2484

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

And of course, the first argument must be x, not f(x), when called with two arguments. So the first call I madethere failed. And if we just want to tell MATLAB the stride alone, again you need to use this form

cumtrapz(1,y)

ans = 1×10

0 4 20 63 154 320 594 1015 1628 2484

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

which works perfectly.

While I know and remember this behavior of cumtrapz and trapz, these tools run counter to our expectations of how the code SHOULD have been written. Well, at least to MY expectations. And given the number of times new and even highly experiencedd users have tripped over this quirk, it is more than just me.

And, yes, I know this is something that could never be fixed, due to the amount of legacy code that relies on tools like trapz and cumtrapz. But it still frustrates me every time I need to use them, because it forces me to waste brain sweat. :)

set(0,'RantMode','off')