18 views (last 30 days)

Show older comments

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

Asvin Kumar
on 12 May 2021

Edited: Asvin Kumar
on 12 May 2021

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

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

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

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

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!