Trapezoidal Rule with symbolic limits
6 views (last 30 days)
Show older comments
Cengizhan Demirbas
on 5 Jan 2021
Answered: Walter Roberson
on 6 Jan 2021
Hello, I am trying to solve this equation by trapezoidal rule from Seppo A. Korpela - Principles of Turbomachinery-Wiley (2011) book:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/479843/image.png)
All variables are known except for K, which I need to find. I defined it as a symbolic variable and tried the following code:
gamma = 4/3; % in the equation
cp = 1148; % in the equation
R = 287;
massflowrate = 3.2; % in the equation
Vx = 150; % in the equation
V1 = Vx;
T1 = 416;
T01 = T1 + (V1^2 / (2*cp)); % in the equation
P1 = 323000;
density1 = P1/(R*T1);
P01 = P1 + (0.5*density1*V1^2);
density01 = P01/(R*T01); % in the equation
alpha2 = 67.18;
V2 = V1/cosd(alpha2);
Vu2 = V2*sind(alpha2); % in the equation
flowcoefficient = 0.5;
U = Vx / flowcoefficient;
shaftspeed = 8320*pi/30;
rm = U/shaftspeed; % in the equation
And the calculation part is:
syms K y
f(y) = (1-Vx^2/(2*cp*T01)-Vu2^2/(2*cp*T01*y^2))^(1/(gamma-1))*y; %only the integral part
a = 2*K/(1+K);
b = 2/(1+K);
n = 10;
h = (b-a)/n;
s = 0.5.*(f(a) + f(b));
for i = 1:n-1
s = s + f(a + i*h);
end
I1 = h*s
eqn1 = I1*2*pi*Vx*rm^2*density01 == massflowrate; %whole equation
hub_to_tip1 = solve(eqn1,K)
Which brings me a 51x1 matrix for hub_to_tip1 with complex roots. My question is that is there a way to solve for K explicitly? The book suggests using numerical integration and that is why I used trapezoidal rule, but perhaps there is a better way to solve this?
4 Comments
Accepted Answer
Walter Roberson
on 6 Jan 2021
syms K positive
syms y
f(y) = (1-Vx^2/(2*cp*T01)-Vu2^2/(2*cp*T01*y^2))^(1/(gamma-1))*y; %only the integral part
a = 2*K/(1+K);
b = 2/(1+K);
I1 = int(f(y),y,a,b);
eqn1 = I1*2*pi*Vx*rm^2*density01 == massflowrate; %whole equation
hub_to_tip1(1) = vpasolve(eqn1, K, 0.1);
hub_to_tip1(2) = vpasolve(eqn1, K, 1);
hub_to_tip1(3) = vpasolve(eqn1, K, 10)
Skip the Simpsons approximation; it is not doing anything useful for you. The int() involves log() so your equation ends up trying to find the root of expressions with variables to power 10 and including log() of the variable. No hope of finding a closed form solution.
Sure, you could taylor it with degree 5, getting out a polynomial of degree 4 and solving for exact roots, but you only get out approximate values and not even enough of them .
0 Comments
More Answers (0)
See Also
Categories
Find more on Calculus 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!