How to remove "root(f(z), z, n)" from solve solution?

Wondering if it is possible to remove the expression "root(f(z), z, n)" from the solution presented by the solve function. In this case, my code is listed below and the outputs of "steadystate = solve(eq5,eq6,eq7,eq8, L,S,B,A)" include the expression "root(f(z), z, n)." I understand it to be the case that this stands in for the nth root of the polynomial f(z), but I'm wondering if it is possible to force a solution from MatLab which does not feature this expression.

clear 

syms f B r L k c g t S u v h m w J a n p V D A mu q m1 m2 m3

assume(f,'positive')
assume(r,'positive')
assume(k,'positive')
assume(c,'positive')
assume(g,'positive')
assume(t,'positive')
assume(u,'positive')
assume(v,'positive')
assume(h,'positive')
assume(m,'positive')
assume(w,'positive')
assume(J,'positive')
assume(a,'positive')
assume(n,'positive')
assume(p,'positive')
assume(V,'positive')
assume(D,'positive')
assume(mu,'positive')
assume(q,'positive')
assume(A,'positive')
assume(m1,'positive')
assume(m2,'positive')
assume(m3,'positive')

eq1 = (f*B + r*L - k*L - c*L - m1*A*L)
eq2 = (g*B + t*S - u*S - v*S - m2*A*S)
eq3 = (h*S + m*L + w*B - a*B - n*B - p*B - m3*A*B)
eq4 = (q*(A+S+B)-mu*A)

eq5 = (eq1 == 0)
eq6 = (eq2 == 0)
eq7 = (eq3 == 0)
eq8 = (eq4 == 0)

steadystate = solve(eq5,eq6,eq7,eq8, L,S,B,A)

[Torsten]

Maybe

steadystate = solve(eq5,eq6,eq7,eq8, L,S,B,A);
steadystate = vpa(steadystate)

[Tristen Jackson]

Unfortunately, using vpa(steadystate) is not effective, nor is applying this function to any of the components of 'steadystate' individually.

Answer 1

[Walter Roberson]

syms f B r L k c g t S u v h m w J a n p V D A mu q m1 m2 m3

assume(f,'positive')
assume(r,'positive')
assume(k,'positive')
assume(c,'positive')
assume(g,'positive')
assume(t,'positive')
assume(u,'positive')
assume(v,'positive')
assume(h,'positive')
assume(m,'positive')
assume(w,'positive')
assume(J,'positive')
assume(a,'positive')
assume(n,'positive')
assume(p,'positive')
assume(V,'positive')
assume(D,'positive')
assume(mu,'positive')
assume(q,'positive')
assume(A,'positive')
assume(m1,'positive')
assume(m2,'positive')
assume(m3,'positive')

eq1 = (f*B + r*L - k*L - c*L - m1*A*L)

eq1 = 

eq2 = (g*B + t*S - u*S - v*S - m2*A*S)

eq2 = 

eq3 = (h*S + m*L + w*B - a*B - n*B - p*B - m3*A*B)

eq3 = 

eq4 = (q*(A+S+B)-mu*A)

eq4 = 

eq5 = (eq1 == 0)

eq5 = 

eq6 = (eq2 == 0)

eq6 = 

eq7 = (eq3 == 0)

eq7 = 

eq8 = (eq4 == 0)

eq8 = 

steadystate = solve(eq5,eq6,eq7,eq8, [L,S,B,A], 'maxdegree', 3)

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

steadystate = struct with fields:
    L: [3×1 sym]
    S: [3×1 sym]
    B: [3×1 sym]
    A: [3×1 sym]

steadystate.L

ans = 

[Walter Roberson]

You might recognize that this is an exact solution, but that for practical purposes it is nearly useless.

If you were planning to use the formulas symbolically in further computation, then you would typically be much better off keeping the root() form, as MATLAB knows how to manipulate those for diff() and int() purposes, and it is a lot faster to do things like take the root() form and cube it than it is to take the exact solution and cube it.

If you were planning to use the formulas numerically, then you would typically be much better off extracting the inside formulas and using sym2poly() and numeric roots()

I have yet to meet anyone who could glance down at a proposed exact root of a cublic and say "That's not right, that + cnv should be a - cnv !" -- the exact formulas are too complex for the great majority of humans to understand.

[Torsten]

Is there any way to "reconstruct" how MATLAB came to the final expressions for the variables ?

[Walter Roberson]

syms A0 A1 A2 A3 x
sols = solve(A3*x^3 + A2*x^2 + A1*x^1 + A0*x^0, x, 'MaxDegree', 3);
string(sols(1))

ans = "(((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3) - (A1/(3*A3) - A2^2/(9*A3^2))/(((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3) - A2/(3*A3)"

string(sols(2))

ans = "(A1/(3*A3) - A2^2/(9*A3^2))/(2*(((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3)) - A2/(3*A3) - (3^(1/2)*((A1/(3*A3) - A2^2/(9*A3^2))/(((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3) + (((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3))*1i)/2 - (((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3)/2"

string(sols(3))

ans = "(3^(1/2)*((A1/(3*A3) - A2^2/(9*A3^2))/(((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3) + (((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3))*1i)/2 - A2/(3*A3) + (A1/(3*A3) - A2^2/(9*A3^2))/(2*(((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3)) - (((A0/(2*A3) + A2^3/(27*A3^3) - (A1*A2)/(6*A3^2))^2 + (A1/(3*A3) - A2^2/(9*A3^2))^3)^(1/2) - A0/(2*A3) - A2^3/(27*A3^3) + (A1*A2)/(6*A3^2))^(1/3)/2"

This shows the pattern; it is the classic solution for roots of a cubic. Everything beyond that is a matter of MATLAB having used coeffs() internally to figure out what A0, A1, A2, A3 match to in the original equation, then substituting.

Well, except that if you are working with polynomials with non-trivial coefficients, it can often be significantly more efficient to ask for the solution to the model polynomial and then subs() in coefficients.

Answer 2

[Torsten]

syms f B r L k c g t S u v h m w J a n p V D A mu q m1 m2 m3
eq1 = (f*B + r*L - k*L - c*L - m1*A*L)
eq2 = (g*B + t*S - u*S - v*S - m2*A*S)
eq3 = (h*S + m*L + w*B - a*B - n*B - p*B - m3*A*B)
eq4 = (q*(A+S+B)-mu*A)
eq5 = (eq1 == 0)
eq6 = (eq2 == 0)
eq7 = (eq3 == 0)
eq8 = (eq4 == 0)
vL = solve(eq5,L)
vS = solve(eq6,S)
vB = subs(eq8,S,vS)
vB = solve(vB,B)
vLL = subs(vL,B,vB)
vSS = subs(vS,B,vB)
vA = subs(eq3,[S,L,B],[vSS,vLL,vB])
simplify(vA)

If you run this code, you will see that you end up with product of poynomials in A of degree 1, 3 and 1.

So the maximum degree is 3 - you can solve analytically for A.

The expressions for B, S and L now follow easily by substituting the expressions for A in vLL, vSS and vB.