How to Solve Non-Linear Equation with changing coefficient

Hello !

Equation at the below is for mach number calculations at different sections in nozzle,

I would like to find mach number

gama = 1.4 
Athroat = 650
(1/mach)*((2*(1+((gama-1)/2)*mach^2))/(gama+1))^((gama+1)/(2*(gama-1)))-Area/Athroat = 0

The thing is that I have to solve this equation for changing Area

Area = [1000 650 800 900 1200] etc.

I tried bisection method, however iteration took like forever.

Is there any short way to solve this equation in Matlab.

Thank you so much.

Answer 1

[Star Strider]

Try this:

gama = 1.4; 
Athroat = 650;
Area = [1000 650 800 900 1200];

machfcn = @(mach,area) (1./mach).*((2*(1+((gama-1)/2).*mach.^2))/(gama+1))^((gama+1)/(2*(gama-1)))-area./Athroat;

for k = 1:numel(Area)
    Mach(k) = fzero(@(mach) machfcn(mach,Area(k)), 10);
end

figure
plot(Area, Mach, 'p')
grid

See the documentation on Anonymous Functions (link)P and fzero (link) for relevant discussions on both.

[Bugra Aksoy]

Thank you so much. this code works. Can I ask Another question ?

Mach(k) = fzero(@(mach) machfcn(mach,Area(k)), 10);

What does 10 imply here ?

Thanks a lot.

[Star Strider]

As always, my pleasure.

The 10 is an initial guess for ‘Mach’. All nonlinear parameter estimation functions need to have an initial estimate for the parameters of interest.

Answer 2

[John D'Errico]

Note there is no need to use fzero. Roots is entirely adequate, and more accurate. Roots also recognizes there are TWO real solutions in general.

syms Area mach
>> gama = 1.4
gama =
          1.4
>> Athroat = 650
Athroat =
   650
>> (1/mach)*((2*(1+((gama-1)/2)*mach^2))/(gama+1))^((gama+1)/(2*(gama-1)))-Area/Athroat
ans =
(mach^2/6 + 5/6)^3/mach - Area/650

It is essentially just a polynomial. I've used syms just to make that clear.

Multiply by mach, fine as long as it is non-zero, and collect coefficients. Multiply by 216 too.

expand(((mach^2/6 + 5/6)^3/mach - Area/650)*mach)
ans =
mach^6 + 15*mach^4 + 75*mach^2 - (108*Area*mach)/325 + 125

So the polynomial coefficients as a function of Area are:

machcoeff = @(Area) [1 0 15 0 75 -108*Area*/325 125];

Now, we can use roots.

roots(machcoeff(1000))
ans =
      0.79428 +     3.8329i
      0.79428 -     3.8329i
      -1.9458 +     2.5675i
      -1.9458 -     2.5675i
       1.8863 +          0i
      0.41673 +          0i

It finds two solutions that are not complex. The others are of no interest to us.

excludecomplex = @(vec) vec(imag(vec) == 0);

Now, put it all together:

excludecomplex(roots(machcoeff(1200)))
ans =
       2.1058
      0.33506

Note there is no solution to the problem when Area is the same as Athroat.

excludecomplex(roots(machcoeff(650)))
ans =
  0×1 empty double column vector

Just loop over the various values for Area.

Only you know which of those two solutions is meaningful, but as long as Area is greater than Athroat, two solutions seem to exist.

Areas =  [660:20:1200];
Machs = NaN(2,length(Areas));
for n = 1:length(Areas)
  M = excludecomplex(roots(machcoeff(Areas(n))));
  if ~isempty(M)
    Machs(:,n) = M;
  end
end
plot(Areas,Machs,'o-')
yline(1);

I imagine one of those solutions is the one you are interested in seeing.

Of course, if you wanted to leave the other coefficients like Athroat and gama in there, you could have done that too.