## How to Solve Non-Linear Equation with changing coefficient

### Bugra Aksoy (view profile)

on 9 Apr 2019
Latest activity Edited by John D'Errico

on 10 Apr 2019

### Star Strider (view profile)

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.

### Star Strider (view profile)

on 9 Apr 2019

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

### Bugra Aksoy (view profile)

on 9 Apr 2019
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

### Star Strider (view profile)

on 9 Apr 2019
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.

### John D'Errico (view profile)

on 10 Apr 2019
Edited by John D'Errico

### John D'Errico (view profile)

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.