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

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...

Answer 1

[John D'Errico]

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.

[Sentient6]

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]

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.