MATLAB Answers

Numerical integration (trapz) between patch

20 views (last 30 days)
Timmy Ngo
Timmy Ngo on 5 May 2021
Edited: Asvin Kumar on 12 May 2021
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

Answers (1)

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

Tags

Community Treasure Hunt

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

Start Hunting!