Numerical solution for a complex formula

5 views (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 (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.
ans = 
Or you could multiply by M, then collect the coefficients of the resulting 6th degree polynominal, then use roots.
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 =
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.





Community Treasure Hunt

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

Start Hunting!