How do fix this code to give me the error?
Show older comments

%% dt=10
h=10;
t=[100:h:1500];
Cp=(.00000000747989.*t.^3)-(.000028903.*t.^2)+.0423484.*t+36.1064;
dH=trapz(t,Cp)
for i=1:5
h=h*10;
tt=100:h:1500;
CcPp=(.00000000747989.*tt.^3)-(.000028903.*tt.^2)+.0423484.*tt+36.1064;
aa=trapz(tt,CcPp);
dHH(i)=aa;
error(i)=abs((dH-dHH(i))/dH);
hh(i)=h;
end
[10 dH 0;hh' dHH' error']
plot(log10(hh),log10(error),'o')
xlabel('Log_{10} (\Delta h)')
ylabel('Log_{10} \epsilon')
a=polyfit(log10(hh),log10(error),1);
f=@(x) a(1)*x+a(2);
hold on
fplot(f,[-5 -1])
hold off
a
%% dt=20
t=[100:20:1500];
Cp=(.00000000747989.*t.^3)-(.000028903.*t.^2)+.0423484.*t+36.1064;
enthalpy_change=trapz(t,Cp)
%% dt=50
t=[100:50:1500];
Cp=(.00000000747989.*t.^3)-(.000028903.*t.^2)+.0423484.*t+36.1064;
enthalpy_change=trapz(t,Cp)
%% dt=100
t=[100:100:1500];
Cp=(.00000000747989.*t.^3)-(.000028903.*t.^2)+.0423484.*t+36.1064;
enthalpy_change=trapz(t,Cp)
%% dt=200
t=[100:200:1500];
Cp=(.00000000747989.*t.^3)-(.000028903.*t.^2)+.0423484.*t+36.1064;
enthalpy_change=trapz(t,Cp)
Answers (2)
Cp = @(T) .00000000747989.*T.^3-.000028903.*T.^2+.0423484.*T+36.1064;
H = @(T) 1/4*.00000000747989.*T.^4-1/3*.000028903.*T.^3+1/2*.0423484.*T.^2+36.1064*T;
deltaT = [10 20 50 100 200];
Tstart = 100;
Tend = 1500;
deltaH_true = H(Tend)-H(Tstart);
for i = 1:numel(deltaT)
dT = deltaT(i);
T = Tstart:dT:Tend;
deltaH(i) = trapz(T,Cp(T));
error(i) = abs((deltaH(i)-deltaH_true)/deltaH_true);
end
format long
error
p = polyfit(log10(deltaT),log10(error),1)
Walter Roberson
on 17 Oct 2022
h=h*10;
tt=100:h:1500;
CcPp=(.00000000747989.*tt.^3)-(.000028903.*tt.^2)+.0423484.*tt+36.1064;
aa=trapz(tt,CcPp);
What is happening is that your increment gets large enough that eventually tt becomes a scalar. But when you pass a scalar to the second parameter of trapz then it is interpreted as a dimension number rather than as a coordinate. (Note that trapz() of a scalar is 0)
Categories
Find more on Numerical Integration and Differentiation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!