How to make sympy only solve the REAL results when solve Higher order equations?

The best way to speed up the solution is probably to change the representation of the problem:

1. make every equations polynomial (applies to f1)
2. make every number in instance of Rational
3. derive a Groebner Basis Representation

I made a jupyter notebook for this, to demonstrate intermediate results, see https://nbviewer.org/github/cknoll/demo-material/blob/main/symbolic/groebner_basis_demo.ipynb

Here is the code (without intermediate results):

from sympy import *
from sympy.solvers.solveset import nonlinsolve

import sympy as sp

def clean_numbers(expr):
    
    numbers = expr.atoms(sp.Number)
    rplmts = [(n, sp.Rational(n)) for n in numbers]
    return expr.subs(rplmts)


h_0 = 170
b = 120
Es = 205
h = 200
M = 19760400
fc = 83
fy = 585


Ec_f, As_f, ft_f = (42.4, 226, 3)
vp_c, vp_s, x_n = xx = symbols('vp_c, vp_s, x_n', real=True)


f1 = vp_c * (h_0 - x_n) - x_n * vp_s # ← converted to a polynomial equation


f2 = 1 / 2 * Ec_f * vp_c * b * x_n - As_f * Es * vp_s - ft_f * b * (h - x_n)
f3 = M - As_f * Es * vp_s * (h_0 - x_n / 3) - ft_f * b * (h - x_n) * (h / 2 + x_n / 6)

ff = clean_numbers(Matrix([f1, f2, f3]).expand())

# reformulate the polynomial equations in terms of a groebner basis
# this set of equations has the same set of solutions but is described by different (often easier) equations
# note that this description depends on the (lexical) order of the variables
gb = groebner(ff, xx, order="lex")

# convert to matrix for easier access
ff2 = Matrix(gb.args[0]) ##:

# solve the last equation (only depends on one variable)

sol2 = solve(ff2[-1], xx[-1]) ##

sol = []
for s in sol2:
    
    tmp_eqns = ff2[:2, :].subs(xx[-1], s.evalf())
    tmp_sol = solve(tmp_eqns, xx[:2]) # → dict
    tmp_sol[xx[-1]] =  s.evalf()
    sol.append(tmp_sol)

print(sol)

→

[{vp_c: 1.00656308634739, vp_s: 2.01548035330581, x_n: 56.6225231688573},
{vp_c: -0.469378305833599, vp_s: 1.02614290514871, x_n: -143.317861965125}]

It takes about 1s.

If you can reduce the system down to a single equation that needs to be solved, you can use real_roots and that should be very fast.

I assume your loop will provide different para_guessed values. Let's just call those x,y,z for now and get the set of equations in terms of those variables. I'm replacing your call to "nonlinsolve" in cal to be return ([f1, f2, f3], [vp_c, vp_s, x_n])

>>> from sympy.abc import x, y, z
>>> eqs, v = cal((x,y,z))

The last two equations are easy to solve for the the first two variables:

>>> v2 = solve(eqs[-2:], v[:2], dict=True)

Now substitute those into the first equation -- there is only one solution so we only need to deal with one possible equation:

>>> xeq = eqs[0].xreplace(v2[0])

To get this most ready to go, let's get the numerator now:

>>> top = xeq.as_numer_denom()[0]

Now... plug in your values -- or any others that you are interested in -- and get solutions for that single variable:

>>> reps = dict(zip((x,y,z),(42.4, 226, 3)))
>>> e = xeq.subs(reps)
>>> xns = real_roots(e)

Backsubstitute to get the other variables:

>>> for i in xns:
...     s = {i:j.n() for i,j in Dict(v2[0]).xreplace(reps).xreplace({v[-1]: i}).items()}
...     s.update({v[-1]: i.n()})
...     print(s)

Gives

{vp_s: 1.02614290514871, vp_c: -0.469378305833599, x_n: -143.317861965125}
{vp_s: 2.01548035330581, vp_c: 1.00656308634739, x_n: 56.6225231688573}

So to recap: you use the CAS to reduce your system symbolically to a single equation for which the real roots can be solve very quickly after substituting in the known values. Then that is used to determine the values of the other equations that were already easy to solve symbolically to get a full numerical solution. No imaginary parts....and quickly.