Not Getting Both Roots to Nonlinear Equation

5 次查看（过去 30 天）
Matthew
Matthew 2017-2-26
评论： Walter Roberson 2017-2-26
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
  2 个评论
Roger Stafford
Roger Stafford 2017-2-26
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
John D'Errico 2017-2-26
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2017-2-26
编辑：John D'Errico 2017-2-26
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

更多回答（2 个）

Walter Roberson
Walter Roberson 2017-2-26
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.
  1 个评论
Walter Roberson
Walter Roberson 2017-2-26
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.

请先登录，再进行评论。

Matthew
Matthew 2017-2-26
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
  1 个评论
Walter Roberson
Walter Roberson 2017-2-26
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by