Not Getting Both Roots to Nonlinear Equation

Can anyone tell me why the following script doesn't give me both answers to the nonlinear equation? The "solve" command only returns one of the equation's roots for M. My TI-89 calculator spits out both roots right away, but I don't know why Matlab is not.

%%Input Section
Ath=0.289;  %Throat area (m^2)
Ae=0.578;   %Exit area (m^2)
gamma=1.37; %Specific heat ratio

%%Calculations
syms M
eqn=Ae/Ath==(1/M)*(((1+((gamma-1)/2)*M^2)/(1+((gamma-1)/2)))^((gamma+1)/2/(gamma-1)));
solM=solve(eqn,M);  %Solve equation for Mach speed
MDbl=double(solM);  %Convert to double numeric

Thanks in advance, M Ridzon

[Roger Stafford]

Consider yourself lucky that you got even one root to that equation. Many equations with fractional powers such as that would receive the message, "Warning: Explicit solution could not be found”.

I would recommend using the numerical solver, ‘fzero’, instead for such equations using different initial guesses.

[John D'Errico]

Roger is correct. Note that when I tried using solve, it took a fairly long time before it gave up and put the problem to vpasolve instead.

Answer 1

[John D'Errico]

No analytical solution exists, so a numerical solver was employed. In general, all roots of a nonlinear equation are not always assured to be found. One can always pose an equation that will mess up any numerical solver.

solM=vpasolve(eqn,M)
solM =
2.1753082190895798956060796154563

If the root is a single root at that location, you can often find the other root by a process called deflation. The idea is to divide by a factor that kills that root off.

vpasolve((((37*M^2)/237 + 200/237)^(237/74)/M - 2)/(M-solM))
ans =
0.30682280437175532207233542049465

Or, if you have plotted the relation, then you would know a decent estimate for the other root.

solM2 = vpasolve(eqn,M,0.2)
solM2 =
0.30682280437175532207233542049465

Answer 2

[Walter Roberson]

The equation is difficult to solve in closed form solution, requiring polynomials of degree 474 and 74th' roots. and fishing around the combinations for the two real solutions. You would be much better off going for vpasolve() one root at a time.

Unfortunately, traditional calculus techniques such as solving for the roots of diff(f(x)^2,x) do not get anywhere useful in this problem -- you can easily find that the maximum is at 1 but finding the 0 crossings symbolically doesn't get far.

[Walter Roberson]

This expression does have solutions for that gamma (and for other gamma expressible as positive rationals). The two solutions are

(1/6309913122)*RootOf(37*Z^474-708843020299236*237^(15/37)*Z^74+141768604059847200*237^(15/37), Z, index = 1)^237*237^(59/74)

(1/6309913122)*RootOf(37*Z^474-708843020299236*237^(15/37)*Z^74+141768604059847200*237^(15/37), Z, index = 2)^237*237^(59/74)

Here, RootOf(F(Z),Z,index=K) stands for the K'th member of the set Z such that F(Z) = 0 . That is, there are 474 distinct roots of the polynomial of degree 474, and the first two roots are the ones you want.

Answer 3

[Matthew]

I guess I'm a little amazed that Matlab struggled with this equation, but my TI-89 calculator didn't hesitate at all. Maybe Texas Instruments and Mathworks need to hook up! Haha!

I really didn't think this would end up being a big deal. Nevertheless, I know what the analytical curve looks like and approximately where the solutions are, so it looks like I'll resort to VPASOLVE.

Thanks, M Ridzon

[Walter Roberson]

When you use a floating point number as an exponent, you abandon most mechanical analytic methods, and do not regain them until you convert the floating point number in to a rational. At that point, though, you may well end up dealing with large powers and large roots. As there is no analytic solution for roots of a polynomial beyond a quartic, the best you can do is get to the well understood "root of a polynomial" representation, such as

RootOf(37*Z^474 - 708843020299236*237^(15/37)*Z^74 + 141768604059847200*237^(15/37), Z)

You can get a general representation along the lines of

RootOf( constant1 * (Z+constant2)^constant3 + constant4, Z )

but unfortunately the constant3 is a fractional power, say p/q. You can expand that out to

RootOf( constant1 * (long expression in Z)^(1/q) + constant4, Z )

You can re-arrange to

constant1 * (long expression in Z)^(1/q) =  -constant4

and raise to the q'th:

constant1^q * (Z+constant2)^(p*q) = (-constant4)^q

or

RootOf( constant1^q * (Z+constant2)^(p*q) - (-constant4)^q , Z)

which can be expressed as a standard polynomial. It will tend to be pretty long, but the terms are predictable by the binomial theorem.

Where does this get you... basically to the point where you could mechanically construct numeric coefficients to pass to root() for numeric evaluation, after which you could filter down to the real-valued solutions.