Simpsons 1/3 Rule to solve integration

13 次查看(过去 30 天)
The enthalpy of a real gas is a function of pressure as described below. The data was taken for a real fluid. Estimate and report the enthalpy of the fluid at 400 K and 50 atm (evaluate the integral from 0.1 atm to 50 atm).
My Code:
% Use Simpsons 1/3 Rule to solve integration
% Lower Limit
a = 0.1; % Pressure in atm
% Upper Limit
b = 50; % Pressure in atm
% Number of Segments
n = 9;
% Declare the function
T = 400; % Temp. in K
V = [250 4.7 2.5 1.49 1.20 0.99 0.75 0.675 0.60]; % Volume at T = 400K
% dVdT = V(T = 350) - V(T = 450)/100
dVdT = [0.625 0.0113 0.005 0.002 0.0014 0.0013 0.001 0.0009 0.0008];
H = V-T.*dVdT;
% inline creates a function of string containing in H
f = inline(H);
% h is the segment size
h = (b - a)/n;
% X stores the summation of first and last segment
X = f(a) + f(b);
% variables Odd and Even to store summation of odd and even terms respectively
Odd = 0;
Even = 0;
for i = 1:2:n - 1
xi = a + (i*h);
Odd = Odd + f(xi);
end
for i = 2:2:n - 2
xi = a+(i*h);
Even = Even + f(xi);
end
% Formula to calculate numerical integration using Simpsons 1/3 Rule
I = (h/3)*(X + 4*Odd + 2*Even)
The Error My Code Produces:
Error using inline (line 51)
Input must be a character vector.
Error in HW6_P1 (line 23)
f = inline(H);
Used this to solve function

回答(2 个)

Cris LaPierre
Cris LaPierre 2022-2-23
As the error states, the input to inline must be a character vector (e.g. f = inline('sin(alpha*x)'))
Note that inline is not recommended. Use Anonymous Functions instead.

Torsten
Torsten 2022-2-23
Be careful here: the pressure values are not equally spaced.
Use
T = 400; % Temp. in K
P = [0.1 5 10 20 25 30 40 45 50]; % Pressure in atm
V = [250 4.7 2.5 1.49 1.20 0.99 0.75 0.675 0.60]; % Volume at T = 400K
% dVdT = V(T = 350) - V(T = 450)/100
dVdT = [0.625 0.0113 0.005 0.002 0.0014 0.0013 0.001 0.0009 0.0008];
H = V-T.*dVdT;
I = trapz(P,H)
to calculate enthalpy at T=400 K and P=50 atm.

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by