Trapezoidal Rule with symbolic limits

Hello, I am trying to solve this equation by trapezoidal rule from Seppo A. Korpela - Principles of Turbomachinery-Wiley (2011) book:

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?

[Walter Roberson]

density1 is not defined at the time it is used.

The next line defines density01 but it is not clear that is the same thing... and it seems unlikely that it is as that would require definining density in terms of itself.

[Cengizhan Demirbas]

Sorry, the original code got kinda messy and i forgot to write the density1, now i edited it in.

[Walter Roberson]

shaftspeed is undefined.

[Cengizhan Demirbas]

Edited it in, hopefully i am not missing any other part

Answer 1

[Walter Roberson]

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 .