Not Getting Both Roots to Nonlinear Equation
1 回表示 (過去 30 日間)
古いコメントを表示
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.
採用された回答
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
0 件のコメント
その他の回答 (2 件)
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
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
2017 年 2 月 26 日
1 件のコメント
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 で Particle Swarm についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!