The answer from Maelstrom is good but I just want to add a few points.
The values you substitute are all floats and with those values the polynomial is ill-conditioned. That means that the form of the expression that you substitute into can affect the accuracy of the returned results. That is one reason why substituting values into the solution from solve does not necessarily give exactly the same value that you get from substituting before calling solve.
Also before you substitute the symbols it isn’t possible for solve to know which of the three roots is real. That’s why you get three solutions from
solve(expr, r) and only one solution from
solve(expr.subs(vdict), r). The third solution which is real after the substitution is the same (ignoring the tiny imaginary part) as returned by solve after the substitution:
In : soln1.subs(vdict).n() Out: 0.000182122477993494 + 1.23259516440783e-32⋅ⅈ In : solve(expr.subs(vdict), r) Out: [0.000182122477993494]
Because the polynomial is ill-conditioned and has a large gradient at the root
nsolve has a hard time finding this root. However
nsolve can find the root if given a narrow enough interval:
In : nsolve(expr.subs(vdict), r, [0.0001821, 0.0001823]) Out: 0.000182122477993494
Since this is essentially a polynomial your best bet is actually to convert it to a polynomial and use
nroots. The quickest way to do this is using
as_numer_denom although in this case that introduces a spurious root at zero:
In : Poly(expr.subs(vdict).as_numer_denom(), r).nroots() Out: [0, 0.000182122477993494, -9.17942953565356e-5 - 0.000158143657514284⋅ⅈ, -9.17942953565356e-5 + 0.000158143657514284⋅ⅈ]
CLICK HERE to find out more related problems solutions.