Message Boards Message Boards

Help solving a cubic equation with a parameter

Posted 3 months ago

POSTED BY: Jackie cao
6 Replies

The third solution is not valid for that value of the parameter:

eq0 = 0 == 
   10125/(14818 \[Phi]) + 
    150 y0 (-6308.297007407409` - 148180.` \[Phi]) + 
    150 (946244.5511111114` + 13000.850814965874`/\[Phi] + 
       1.11135`*^7 \[Phi] - (
       0.022395039313571143` (6959.` + (337500.` - 
             2250.` y0) \[Phi])^(3/2))/\[Phi]);
\[Phi]0 = 0.00145;
Solve[eq0, y0, Assumptions -> \[Phi] == \[Phi]0]
Solve[eq0 /. \[Phi] -> \[Phi]0, y0]
sol3 = Solve[eq0, y0] /. \[Phi] -> \[Phi]0
eq0[[2]] /. \[Phi] -> \[Phi]0 /. sol3

The formulas for the first and second solution contain an imaginary unit, but it is multiplied by a very small coefficient, probably negligible. The main problem is that all three solutions have fractional powers, which are problematic.

POSTED BY: Gianluca Gorni

Yes, it can be ignored. This is the casus irreducilbilis, in which the solutions to a cubic can be expressed only in terms of n-th roots that have complex values. The small imaginary parts are the result of round-off error.

If you'd like to get rid of them and don't mind Root[] objects, specify Cubics -> False on the exact system:

Solve[
   Rationalize[Rationalize[
     0 == 
      10125/(14818 \[Phi]) + 
       150 y0 (-6308.297007407409` - 148180.` \[Phi]) + 
       150 (946244.5511111114` + 13000.850814965874`/\[Phi] + 
          1.11135`*^7 \[Phi] - (
          0.022395039313571143` (6959.` + (337500.` - 
                2250.` y0) \[Phi])^(3/2))/\[Phi])
     ], 0]
   , {y0}, Cubics -> False
   ] /. Rationalize[\[Phi] -> 0.00145] // N

Solve::nongen: There may be values of the parameters for which some or all solutions are not valid.

(* {{y0 -> -13.9237}, {y0 -> 13.8775}, {y0 -> 1712.28}}  *)

Note that in this case Solve[] warns us that not all solutions may be valid, which indeed turned out to be the case. Root[] computes the real roots of a polynomial as Real numbers without complex round-off errors.

If you want a more rigorous method, use Method -> Reduce and increase MaxExtraConditions. Here MaxExtraConditions -> 1 is enough, but I started with MaxExtraConditions -> Infinity:

solRat = Solve[
   Rationalize[Rationalize[
     0 == 
      10125/(14818 \[Phi]) + 
       150 y0 (-6308.297007407409` - 148180.` \[Phi]) + 
       150 (946244.5511111114` + 13000.850814965874`/\[Phi] + 
          1.11135`*^7 \[Phi] - (
          0.022395039313571143` (6959.` + (337500.` - 
                2250.` y0) \[Phi])^(3/2))/\[Phi])
     ], 0]
   , {y0}, Cubics -> False, Method -> Reduce, MaxExtraConditions -> 1
   ]; (* 7 solutions depending on cases for the parameter phi *)
solPhi =(* select the defined cases *)
  Select[% /. Rationalize[\[Phi] -> 0.00145] // RootReduce, 
   NumericQ[y0 /. #] &];
solPhi // N (* we see the original was correct *)

Solve::useq: The answer found by Solve contains equational condition(s) {0==(192677381803147327908319230731396042083840204210666350+<<22>>)/(246032254379437568608800 (9387804582563199208933776766+<<1>>+138587276555703258030203175 Root[<<1>>&,1])),0==<<1>>/(2<<21>>00<<1>><<1>>),0==<<1>>/(246032254379437568608800<<1>><<1>>)}. A likely reason for this is that the solution set depends on branch cuts of Wolfram Language functions.

(*  {{y0 -> -13.9237}, {y0 -> 13.8775}}  *)

This time we get more of an idea why the previous Solve[] said not all solutions may be valid. We also see that several other solution-cases were rejected by the OP's original Solve[].

In the OP's second code, increasing MaxExtraConditions yields conditions on the solutions that are helpful in filtering out bogus solutions:

solApprox = 
  Solve[0 == 
    10125/(14818 \[Phi]) + 
     150 y0 (-6308.297007407409` - 148180.` \[Phi]) + 
     150 (946244.5511111114` + 13000.850814965874`/\[Phi] + 
        1.11135`*^7 \[Phi] - (
        0.022395039313571143` (6959.` + (337500.` - 
              2250.` y0) \[Phi])^(3/2))/\[Phi]), y0,
   Cubics -> False, MaxExtraConditions -> 1];
solApprox /. \[Phi] -> 0.00145
(* <warnings omitted>
{{y0 -> Undefined}, {y0 -> Undefined}, {y0 -> Undefined}, {y0 -> Undefined},
 {y0 -> -13.9237}, {y0 -> 13.8775}, {y0 -> Undefined}}
*)

If we remove the conditions from the solutions with Normal[], we see the last three are the solutions returned by the OP's second Solve[]:

Normal@solApprox /. \[Phi] -> 0.00145
(*
{{y0 -> -68.88}, {y0 -> -67.7773}, {y0 -> 41.1304 - 66.9167 I}, {y0 -> 41.1304 + 66.9167 I},
 {y0 -> -13.9237}, {y0 -> 13.8775}, {y0 ->1712.28}}
*)
POSTED BY: Michael Rogers

Thank you for your help. Yeah, you're right. I realized the real part of the first and second solution is the right answer. In the past two days I have scrutinized the general solution of cubic equations with parameters in search of why they carry imaginary units, and have discovered a possible range of values problem and a precision problem. The precision was too high resulting in imaginary numbers being left in the equation that would have cancelled out. In addition, the 1.5th power in the equation also caused the third solution to carry imginary units, and this is what gave me a headache. But I can't change that this power of a fraction is also integrated over from a power of a fraction. /facepalm

POSTED BY: Jackie cao

But I wonder if even a small imaginary unit can be ignored?

POSTED BY: Jackie cao

When an equation has floating-point coefficients, the concept of solution is a little fuzzy. Try replacing into the original equation, with and without the small imaginary component. Can you see much difference?

POSTED BY: Gianluca Gorni
Posted 2 months ago

Thank you for these Instructions and symbols, this will be helpful to solve polynomial parametric equations.

POSTED BY: Updating Name
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract