Not Getting Both Roots to Nonlinear Equation

Question

Matthew 2017-2-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/326905-not-getting-both-roots-to-nonlinear-equation

评论： 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 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 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.

请先登录，再进行评论。

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

Answer 1

John D'Errico 2017-2-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/326905-not-getting-both-roots-to-nonlinear-equation#answer_256317

编辑：John D'Errico 2017-2-26

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Walter Roberson 2017-2-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/326905-not-getting-both-roots-to-nonlinear-equation#answer_256315

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2017-2-26

在 MATLAB Online 中打开

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 2017-2-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/326905-not-getting-both-roots-to-nonlinear-equation#answer_256346

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2017-2-26

在 MATLAB Online 中打开

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.

请先登录，再进行评论。