"Solve" function and "overly precise" parameters??

1 次查看（过去 30 天）
Sentient6
Sentient6 2018-4-6
评论： John D'Errico 2018-4-7
So I've encountered a really weird problem trying to solve an implicit equation:
eq = eps == (1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1)));
M1 = double(solve(eqi,Mi));
"eps" and "g" are given parameters, and Mi is what I'm trying to calculate. The thing is, when I use "g" with a single digit after a decimal point (I started with 1.2, but later found out 1.1, 1.3, 1, 2, etc), the code works like a charm. But when I try to be more precise (g = 1.21, 1.19, 1.2001), Matalb fails to solve the equation and gives me a strange value that is definitely incorrect. I also get a "warning: Cannot solve symbolically. Returning a numeric approximation instead. " orange error (it doesn't stop running).
It's the first time I've encountered something so weird, I can't even start to guess what the problem might be...

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

采纳的回答

John D'Errico
John D'Errico 2018-4-6
编辑：John D'Errico 2018-4-6
Look more carefully at what results. When g is fairly simple, i.e.,
g = 1.1;
pretty((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1))))
/ 2 \21/2
| Mi 20 |
| --- + -- |
\ 21 21 /
----------------
Mi
We can resolve that by squaring it. Even better, if g is an integer,
g = 1;
pretty((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1))))
1
--
Mi
So when g is pretty basic, the resulting polynomial in Mi is also simple.
But look at what happens when g has a few more digits beyond the decimal point in it.
g = 1.273;
pretty((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1))))
/ 2 \2273/546
| 273 Mi 2000 |
| ------- + ---- |
\ 2273 2273 /
--------------------------
Mi
You MIGHT decide to essentially raise that expression to the 546'th power, the result being a polynomial in Mi of order 2*2273=4546.
  2 个评论
Sentient6
Sentient6 2018-4-7
Well damn, I see it now. Thanks for the answer.
A follow-up - is there a quick alternative way to solve it for more complex 'g' values? A different function than "solve"?
John D'Errico
John D'Errico 2018-4-7
Hmm. My answer was easy. No thought involved. Now you want me to think? ;-)
First, there appear to be two positive solutions for any value of eps, except 1, where there will be a double root, at Mi = 1.
subs((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1))),Mi,1)
ans =
1
The min happens at Mi == 1, for any g.
g seems to be a number greater than 1 in your question, so I''ll plot things for g=1.3. And, since things seem to get pretty large quite fast, I'll implicitly get rid of that exponent, by taking the log.
g = 1.3;
ezplot(log((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1)))),[0 5])
grid on
By the way, using eps as a variable is just a flat out bad idea, given the value of the function eps in MATLAB.
So it looks like things are pretty tame now. Depending on the constant you have for log(eps), and which root you are interested in, this should be easily solvable using numerical methods. However, I hardly expect you will ever have an analytical solution for any general g.
If I wanted to solve this using numerical methods, I can think of a few ways I might go. fzero or vpasolve would be the obvious.
E = 3;
vpasolve(log((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1)))) - log(E),2)
ans =
2.5140963815086844614466043240907
vpasolve(log((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1)))) - log(E),.5)
ans =
0.19958189654964945539477389458526
To be honest, for either root, you could construct a simple fixed point iteration loop that would converge quite nicely, but why bother, since fzero will be easier.
fzero(@(Mi) log((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1)))) - log(E),.5)
ans =
0.19958189654965
fzero(@(Mi) log((1/Mi)*((2/(g+1))*(1+ ((g-1)/2)*Mi^2))^((g+1)/(2*(g-1)))) - log(E),2)
ans =
2.51409638150868
I doubt you will find an analytical solution though.

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by