Numerical integration (trapz) between patch

Hello! I want to numerically integrate two peaks using the trapezoidal method, but I'm afraid that by using the trapz function I'm not integrating over the red marked area (see figure). My suspicion is that the first integrated peak should have been positive. Is there some way around this?

The attached file is required to run the code, I imported the .txt file and Output Type is set to Colum vectors:

close
temp = sgolayfilt(temp,1,51); %Savitzky-Golay fit
H = sgolayfilt(H,1,51); %Savitzky-Golay fit
temp1 = temp(1050:1539); H1 = H(1050:1539); %Interval 1st peak
temp2 = temp(2450:3050); H2 = H(2450:3050); %Interval 2nd peak
I1 = trapz(temp1,H1) %Trapezoidal integration 1st peak
I2 = trapz(temp2,H2) %Trapezoidal integration 2nd peak

figure('Position',[500 500 800 490])
plot(temp,H,'k','LineWidth',1.1)
grid minor
hold on
% area(temp1,H1,'FaceColor','r') "trapz(temp1,H1)" integrates over this?
% area(temp2,H2,'FaceColor','r') "trapz(temp2,H2)" integrates over this?
patch(temp1,H1,'red')
patch(temp2,H2,'red')
xlabel('$T$  [C$^\circ$]','interpreter','latex')
ylabel('$\Delta H$  [mW]','interpreter','latex')
title('PET -- Differential scanning calorimetry curve, where $T\in[0,350^\circ$ C].' ,'interpreter','latex')
legend({'DSC Curve','1st Peak','2nd Peak'},'Location','southwest','interpreter','latex')
xlim([0,360])

I1 = -71.3474

I2 = -156.3513

Many thanks,

Timmy

Answer 1

[Asvin Kumar]

The reason I1 and I2 are negative is because the H values are negative. trapz numerically integrates over the specified range by approximating segments as trapeziums. The heights of these trapeziums is measured from y=0. Since H is negative, all your trapeziums have negative height. Click the hyperlink to learn more about trapz calculation.

It's easy to fix the computation. You just have to subtract the extra area that was included in trapz's computation. Have a look at the script below which shows what area you expected it to integrate vs what it actually integrated over.

Here's a modified version of your script. The variables Peak1 and Peak2 should contain the areas that you're looking for.

%% Set up the Import Options and import the data
opts = delimitedTextImportOptions("NumVariables", 3);
% Specify range and delimiter
opts.DataLines = [3, 4047];
opts.Delimiter = "\t";
% Specify column names and types
opts.VariableNames = ["t", "temp", "H"];
opts.VariableTypes = ["double", "double", "double"];
% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Import the data
DSCPET = readtable("DSCPET.txt", opts);
temp = DSCPET.('temp');
H = DSCPET.('H');
clear opts
%% User's Script
temp = sgolayfilt(temp,1,51); %Savitzky-Golay fit
H = sgolayfilt(H,1,51); %Savitzky-Golay fit
temp1 = temp(1050:1539); H1 = H(1050:1539); %Interval 1st peak
temp2 = temp(2450:3050); H2 = H(2450:3050); %Interval 2nd peak
I1 = trapz(temp1,H1) %Trapezoidal integration 1st peak

I1 = -71.3099

I2 = trapz(temp2,H2) %Trapezoidal integration 2nd peak

I2 = -156.3930

%% User's Plot
% figure('Position',[500 500 800 490])
plot(temp,H,'k','LineWidth',1.1)
grid minor
hold on
% area(temp1,H1,'FaceColor','r') "trapz(temp1,H1)" integrates over this?
% area(temp2,H2,'FaceColor','r') "trapz(temp2,H2)" integrates over this?
patch(temp1,H1,'red')
patch(temp2,H2,'red')
xlabel('$T$  [C$^\circ$]','interpreter','latex')
ylabel('$\Delta H$  [mW]','interpreter','latex')
title('PET -- Differential scanning calorimetry curve, where $T\in[0,350^\circ$ C].'...
    ,'interpreter','latex')
legend({'DSC Curve','1st Peak','2nd Peak'},'Location','southwest',...
    'interpreter','latex')
xlim([0,360])

%% Modified Calculation - Subtract extra trapezium area whose base is at y=0
Peak1 = I1 - 0.5*(H1(1)+H1(end))*(temp1(end)-temp1(1)) 

Peak1 = 19.9070

Peak2 = I2 - 0.5*(H2(1)+H2(end))*(temp2(end)-temp2(1))

Peak2 = -35.4004

%% Modified plot - This is what trapz calculates
figure
plot(temp,H,'k','LineWidth',1.1)
grid minor
hold on
patch([temp1(1); temp1; temp1(end)],[0; H1; 0],'g')
patch([temp2(1); temp2; temp2(end)],[0; H2; 0],'g')
xlabel('$T$  [C$^\circ$]','interpreter','latex')
ylabel('$\Delta H$  [mW]','interpreter','latex')
title('PET -- Differential scanning calorimetry curve, where $T\in[0,350^\circ$ C].'...
    ,'interpreter','latex')
legend({'DSC Curve','Trapz 1','Trapz 2'},'Location','southwest',...
    'interpreter','latex')
xlim([0,360])