MATLAB Answers

0

How to Solve Non-Linear Equation with changing coefficient

Asked by Bugra Aksoy on 9 Apr 2019
Latest activity Edited by John D'Errico
on 10 Apr 2019
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.

  0 Comments

Sign in to comment.

2 Answers

Answer by Star Strider
on 9 Apr 2019
 Accepted Answer

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.

  2 Comments

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.
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.

Sign in to comment.


Answer by John D'Errico
on 10 Apr 2019
Edited by John D'Errico
on 10 Apr 2019

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.

  0 Comments

Sign in to comment.