MATLAB Answers

Numerical integration (trapz) between patch

조회 수: 7(최근 30일)
Timmy Ngo
Timmy Ngo 2021년 5월 5일
편집: Asvin Kumar 2021년 5월 12일
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

답변(1개)

Asvin Kumar
Asvin Kumar 2021년 5월 12일
편집: Asvin Kumar 2021년 5월 12일
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])

태그

Community Treasure Hunt

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

Start Hunting!

Translated by