Trying to solve for P in the equation. Been told fzero does not give correct solution.

%clc; clear; close all; format short g; 
%knowns
gamma = 1.66;
R = 8.314; 
m = 4.0e-3;
RoM = R/m;
P0 = 10e5; 
T0 = 300;
r0 = 0.01;
L = 0.05;
rho0 = P0/(RoM*T0);
 
%function for area
x = [0:0.001:2*L];
func = @(x) r0*((2*x./L)+ exp(-2*x./L));
y = func(x);
plot(x,y);
At = pi*r0*r0;

% equations
mc = At*(sqrt(gamma*P0*rho0))*(2/(gamma+1))^((gamma+1)/(2*(gamma-1))); %mass flux choking
 
%solving for P
f = @(P) ((P/P0)^(2/gamma)) - ((P/P0)^(1+(1/gamma))) == ((gamma-1)/(2*gamma))*((mc^2)/((At^2)*P0*rho0));
r = solve (f,eps)

[dpb]

Is there a specific Q? here?

Have you plotted the function to investigate behavoir???

Answer 1

[Samatha Aleti]

Hi,

I think you need to specify the variables as “sym” to use “solve”. You may also solve for “P” using “fsolve” function.

Here is the documentation link for “fsolve”.

[Walter Roberson]

fsolve cannot deal with the == so you will need to take one side minus the other side and fsolve the resulting expression.

[dpb]

Yeah, Walter. And I tried messing around w/ it some but couldn't get a set of numbers to actually cross zero...it gets pretty close at a number only about half the P0 value...

You first have to vectorize and rewrite the functional to remove the "==" as you note:

>> fnf
fnf =
  function_handle with value:
    @(P)((P./P0).^(2./gamma))-((P./P0).^(1+(1./gamma)))-((gamma-1)./(2.*gamma)).*((mc.^2)./((At.^2).*P0.*rho0))
>> 

Then

P=linspace(P0/10,P0*10,1000);
e=fnf(P);
plot(P,e)
hAx=gca;
hAx.XScale='log'
ylim([-0.5 0.1])

results in

>> [emx imx]=max(e);
>> P(imx), P(imx)/P0
ans =
   4.8649e+05
ans =
    0.4865
>> 

No idea whether this is what OP expects should look like or not; never came back...

fsolve without more work does fail because there isn't a zero crossing and the slope in the region is small--

>> fsolve(fnf,P0/2)

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>
ans =
      500000
>> 

Same thing happens at every starting point I tried but I didn't take time to mess with internals.

[Walter Roberson]

If you convert everything to rational, except using the transcendentatal pi instead of an approximation of it, then 281250*sqrt(3) becomes an exact solution. That is approximately 487139.2898.

[dpb]

Ah, indeed! I hadn't really looked at the algebra, just the numerical solution. Nicely done, Walter!