Numerical solution for a complex formula

1 view (last 30 days)
Owen McKim
Owen McKim on 24 Nov 2021
Commented: John D'Errico on 24 Nov 2021
I am trying to rearrange a complex formula to give a numerical output.
A = ((1)/(M))*((1+((y-1)/(2))*(M^(2)))/((y+1)/(2)))^((y+1)/(2*(y-1)))
This is the formula I am trying to rearrange, to rearrange for M, using the solve function.
Fortunately, the Y is a constant value of 1.4. Therefore, I simplified the formula in MATLAB.
solve((A) == (1/M)*(((1+((1.4-1)/2)*(M^2))/((1.4+1)/2))^((1.4+1)/(2*(1.4-1)))),A)
A = ((M^2/6 + 5/6)^3)/M
I then tried using the solve function but solving for M with the simplied formula but only got imaginary roots returned. I then tried subbing A as a value, 1.2, and still, MATLAB returned the imaginary root.
Now using 1.2 = ((M^2/6 + 5/6)^3)/M to solve for M. I tried using symbolab.com (An online calculator) And that gave me 2 answers of M = 0.59024 and M= 1.53414. These values are correct, as I confirmed from a table.
So how can I get MATLAB to output these 2 M values and not roots? From 1.2 = ((M^2/6 + 5/6)^3)/M

Answers (1)

John D'Errico
John D'Errico on 24 Nov 2021
syms M A
y = 1.4;
eqn = A == ((1)/(M))*((1+((y-1)/(2))*(M^(2)))/((y+1)/(2)))^((y+1)/(2*(y-1)))
eqn = 
Now you wish to solve that expression for M. Note that this expression is equivalent to a 6th degree polynomial, As such, it most likely cannot be resolved into a set of algebraically represented roots. So until you have a value for A, then you can do essentially nothing. (Occasionally you can find a polynomial of degree higher than 4 with algebraicly representable roots, biut that is the exception rather than the rule.)
But if A has a numerical value, then This is a 6th degree polynonial with constant coefficients. So you can use tools to solve it numerically.
vpasolve(subs(eqn,A,1.2))
ans = 
Or you could multiply by M, then collect the coefficients of the resulting 6th degree polynominal, then use roots.
  2 Comments
John D'Errico
John D'Errico on 24 Nov 2021
Just delete the elements of that result that have a non-zero imaginary part.
result = vpasolve(subs(eqn,A,1.2));
result(imag(result) ~= 0) = []
result =
0.59024876098845097803135562613058
1.5341497672017611068224700524926
Often safer when working with floating point numbers will be to compare the imaginary part to a tiny number to decide if you should delete it.
result(abs(imag(result)) < eps) = []

Sign in to comment.

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!