Solve of cubic and quartic equations too slow

Question

In my problem I have a third order algebraic equation for the variable sigma, all other letters are parameters. Here is it's right-hand side, left hand side is zero (you can copy it into your programm):

eq = Er k^2 γ ξ (1 + λ Cos[2 ϕ]) (3 I k Cos[ϕ] + 4 σ Cos[2 ϕ] + 
    I k Cos[3 ϕ]) Subscript[α, ac] + (ξ + σ) (σ + I k Cos[ϕ]) (8 k^4 + 
    16 k^2 tf1 + 2 k^4 γ + k^4 γ λ^2 + 8 Er k^2 γ σ + 16 Er tf1 γ σ + 
    4 k^4 γ λ Cos[2 ϕ] + k^4 γ λ^2 Cos[4 ϕ] + 
    4 k^2 (k^2 + Er γ σ) Sin[2 ϕ]^2 Subscript[μ, 1] + 
    4 (k^2 + Er γ σ) (k^2 + tf1 - tf1 Cos[2 ϕ]) Subscript[μ, 2])

My purpose here is to get Taylor expansion in parameter k for all the roots near k=0:

sol = Solve[eq == 0, σ]
sol = σ /. sol
Series[sol[[1]], {k, 0, 2}

The thing is, that evalution of the line sol = Solve[eq == 0, σ] takes really long time. It takes couple minutes here, but when I am trying to solve similar equation of the 4th order, it takes forever. Note, that in the line sol = Solve[eq == 0, σ], I didn't ask to simplify. But looks like mathematica silently simplifies it, since it is running for a long time. Note, that for the third and forth order equations exact algebraic expressions for the solutions are availible, and indeed, if I type something like Solve[a*x^3+b*x^2+c*x+d==0,x], it gives me the answer instantly. It should be really fast procedure to substitute corresponding expressions into the final formula.

Is there a way to get Taylor expansion I need faster? Solving the equations takes much longer time, than Taylor series.

In case someone needs, this is the 4th order equation, which I cannot solve even in 30 minutes:

-Er k^2 γ ξ (1 + λ Cos[2 ϕ]) (4 (k^2 + 2 σ (σ + τ1)) Cos[2 ϕ] + k^2 (3 + Cos[4 ϕ]))
    Subscript[α, ac] - (ξ + σ) (k^2 + 2 σ (σ + τ1) + k^2 Cos[2 ϕ]) (8 k^4 + 16 k^2 tf1 + 
    2 k^4 γ + k^4 γ λ^2 + 8 Er k^2 γ σ + 16 Er tf1 γ σ + 4 k^4 γ λ Cos[2 ϕ] + 
    k^4 γ λ^2 Cos[4 ϕ] + 4 k^2 (k^2 + Er γ σ) Sin[2 ϕ]^2 Subscript[μ, 1] + 
    4 (k^2 + Er γ σ) (k^2 + tf1 - tf1 Cos[2 ϕ]) Subscript[μ, 2])

Evaluation of sol = Solve[eq == 0,σ] takes forever.

I found out, that Defer can be used to supress evalution, but I need a different thing.

Thanks, MIkhail

Answer 1

Here is a workaround:

eq = Er k^2 γ ξ (1 + λ Cos[2 ϕ]) (3 I k Cos[ϕ] + 4 σ Cos[2 ϕ] + 
    I k Cos[3 ϕ]) Subscript[α, ac] + (ξ + σ) (σ + I k Cos[ϕ]) (8 k^4 + 
    16 k^2 tf1 + 2 k^4 γ + k^4 γ λ^2 + 8 Er k^2 γ σ + 16 Er tf1 γ σ + 
    4 k^4 γ λ Cos[2 ϕ] + k^4 γ λ^2 Cos[4 ϕ] + 
    4 k^2 (k^2 + Er γ σ) Sin[2 ϕ]^2 Subscript[μ, 1] + 
    4 (k^2 + Er γ σ) (k^2 + tf1 - tf1 Cos[2 ϕ]) Subscript[μ, 2]);
coeffs = Reverse@CoefficientList[eq, σ];
sol = σ /. Solve[FromDigits[Array[#[] &, Length[coeffs], 1],
 σ] == 0, σ] /. n_Integer[] :> coeffs[[n]];
Series[sol[[1]], {k, 0, 2}]

Instead of putting eq in Solve, I put a polynomial with dummy variables (1[], 2[], etc.) and then replaced the dummies with the corresponding coefficients in eq.

It seems that this method takes about 3 seconds on my computer when I first ran it; 0.03 seconds for subsequent iterations.

The bottom 3 lines (from coeffs = to the end of code) should work for any equation eq.

Answer 2

I would go about this differently. You have, in effect, a curve parametrized by sigma and k (regard the other parameters as fixed for the moment). This curve is the zero set for eq. You are interested in expanding sigma as a series in k, around k=0. This can be done by setting up equations for derivatives of sigma regarded explicitly as a function of k; these equations arise from differentiating the original one. As we are evaluating derivatives at k=0 we can make that substitution after differentiating but before solving. Here is the process.

eq = Er k^2 γ ξ (1 + λ Cos[
        2 ϕ]) (3 I k Cos[ϕ] + 4 σ Cos[2 ϕ] + 
      I k Cos[3 ϕ]) Subscript[α, 
     ac] + (ξ + σ) (σ + I k Cos[ϕ]) (8 k^4 + 
      16 k^2 tf1 + 2 k^4 γ + k^4 γ λ^2 + 
      8 Er k^2 γ σ + 16 Er tf1 γ σ + 
      4 k^4 γ λ Cos[2 ϕ] + 
      k^4 γ λ^2 Cos[4 ϕ] + 
      4 k^2 (k^2 + Er γ σ) Sin[
         2 ϕ]^2 Subscript[μ, 1] + 
      4 (k^2 + Er γ σ) (k^2 + tf1 - 
         tf1 Cos[2 ϕ]) Subscript[μ, 2]);
eqk = eq /. σ -> σ[k];
derivsys = Table[D[eqk, {k, j}], {j, 0, 3}];
derivvars = Table[D[σ[k], {k, j}], {j, 0, 3}];
newvars = Array[d, 4, 0];
subs = Thread[derivvars -> newvars];
sys = (derivsys /. subs) /. k -> 0;
soln = Solve[sys == 0, newvars];

During evaluation of In[425]:= Solve::svars: Equations may not give solutions for all "solve" variables. >>

In[436]:= Most[newvars].(k^Range[0, 2]/Factorial[Range[0, 2]]) /. soln

(* Out[436]= {(k^2 Sec[ϕ] (16 tf1 Cos[ϕ] + 
      3 Er γ Cos[ϕ] Subscript[α, ac] + 
      3 Er γ λ Cos[ϕ] Cos[
        2 ϕ] Subscript[α, ac] + 
      Er γ Cos[3 ϕ] Subscript[α, ac] + 
      Er γ λ Cos[2 ϕ] Cos[
        3 ϕ] Subscript[α, ac] + 
      4 tf1 Cos[ϕ] Subscript[μ, 2] - 
      4 tf1 Cos[ϕ] Cos[2 ϕ] Subscript[μ, 
       2]))/(4 Er tf1 γ (-4 - Subscript[μ, 2] + 
      Cos[2 ϕ] Subscript[μ, 2])), -I k Cos[ϕ] + (
  k^2 (1 + λ Cos[2 ϕ]) (-3 Cos[ϕ] + 
     4 Cos[ϕ] Cos[2 ϕ] - 
     Cos[3 ϕ]) Sec[ϕ] Subscript[α, ac])/(
  4 tf1 (-4 - Subscript[μ, 2] + 
     Cos[2 ϕ] Subscript[μ, 2])), -ξ - (
  k^2 Cos[2 ϕ] (1 + λ Cos[2 ϕ]) Subscript[α, 
   ac])/(tf1 (-4 - Subscript[μ, 2] + 
     Cos[2 ϕ] Subscript[μ, 2]))} *)

Astute readers will notice there was a wrinkle. Due to multiplicity and/or factors of sigma in the reduced derivative equations (by reduced I mean with the k->0 substitution), we get underdetermined systems. A way around this is to take more derivatives than we actually need, enough to get solutions for the lower derivatives that we require in the expansion. Above I used one extra derivative. I'm sure there is a theory behind this need for extra derivatives (some aspect of prolongation, maybe) but I do not know the details at all.

Answer 3

You can use the new in M12 function AsymptoticSolve for this:$\unicode{958}$

eq = Er k^2 γ ξ (1 + λ Cos[2 ϕ]) (3 I k Cos[ϕ] + 4 σ Cos[2 ϕ] + 
I k Cos[3 ϕ]) Subscript[α, ac] + (ξ + σ) (σ + I k Cos[ϕ]) (8 k^4 + 
16 k^2 tf1 + 2 k^4 γ + k^4 γ λ^2 + 8 Er k^2 γ σ + 16 Er tf1 γ σ + 
4 k^4 γ λ Cos[2 ϕ] + k^4 γ λ^2 Cos[4 ϕ] + 
4 k^2 (k^2 + Er γ σ) Sin[2 ϕ]^2 Subscript[μ, 1] + 
4 (k^2 + Er γ σ) (k^2 + tf1 - tf1 Cos[2 ϕ]) Subscript[μ, 2]);

AsymptoticSolve[eq == 0, σ, {k, 0, 2}] //TeXForm

$\left\{\left\{\sigma \to \frac{k^2 \left(3 \gamma \operatorname{Er} \lambda \alpha _{\operatorname{ac}} \cos (2 \phi )+\gamma \operatorname{Er} \lambda \alpha _{\operatorname{ac}} \cos (2 \phi ) \cos (3 \phi ) \sec (\phi )+\gamma \operatorname{Er} \alpha _{\operatorname{ac}} \cos (3 \phi ) \sec (\phi )+3 \gamma \operatorname{Er} \alpha _{\operatorname{ac}}-4 \mu _2 \operatorname{tf1} \cos (2 \phi )+4 \mu _2 \operatorname{tf1}+16 \operatorname{tf1}\right)}{4 \gamma \operatorname{Er} \operatorname{tf1} \left(\mu _2 \cos (2 \phi )-\mu _2-4\right)}\right\},\left\{\sigma \to \frac{k^2 \left(4 \lambda \alpha _{\operatorname{ac}} \cos ^2(2 \phi )-3 \lambda \alpha _{\operatorname{ac}} \cos (2 \phi )-\lambda \alpha _{\operatorname{ac}} \cos (3 \phi ) \cos (2 \phi ) \sec (\phi )+4 \alpha _{\operatorname{ac}} \cos (2 \phi )-\alpha _{\operatorname{ac}} \cos (3 \phi ) \sec (\phi )-3 \alpha _{\operatorname{ac}}\right)}{4 \operatorname{tf1} \left(\mu _2 \cos (2 \phi )-\mu _2-4\right)}-i k \cos (\phi )\right\},\left\{\sigma \to \frac{k^2 \left(-\lambda \alpha _{\operatorname{ac}} \cos ^2(2 \phi )-\alpha _{\operatorname{ac}} \cos (2 \phi )\right)}{\operatorname{tf1} \left(\mu _2 \cos (2 \phi )-\mu _2-4\right)}-\xi \right\}\right\}$

I tested this result against JungHwan's result for some random values of the parameters, and they agree. The above output is much more concise than JungHwan's results.

Update

Until you have M12, you can use CloudEvaluate:

$Version
CloudEvaluate[System`AsymptoticSolve[eq==0, σ, {k, 0, 2}]] //TeXForm

"11.3.0 for Mac OS X x86 (64-bit) (March 6, 2018)"

$\left\{\left\{\sigma \to \frac{k^2 \left(3 \gamma \operatorname{Er} \lambda \alpha _{\operatorname{ac}} \cos (2 \phi )+\gamma \operatorname{Er} \lambda \alpha _{\operatorname{ac}} \cos (2 \phi ) \cos (3 \phi ) \sec (\phi )+\gamma \operatorname{Er} \alpha _{\operatorname{ac}} \cos (3 \phi ) \sec (\phi )+3 \gamma \operatorname{Er} \alpha _{\operatorname{ac}}-4 \mu _2 \operatorname{tf1} \cos (2 \phi )+4 \mu _2 \operatorname{tf1}+16 \operatorname{tf1}\right)}{4 \gamma \operatorname{Er} \operatorname{tf1} \left(\mu _2 \cos (2 \phi )-\mu _2-4\right)}\right\},\left\{\sigma \to \frac{k^2 \left(4 \lambda \alpha _{\operatorname{ac}} \cos ^2(2 \phi )-3 \lambda \alpha _{\operatorname{ac}} \cos (2 \phi )-\lambda \alpha _{\operatorname{ac}} \cos (3 \phi ) \cos (2 \phi ) \sec (\phi )+4 \alpha _{\operatorname{ac}} \cos (2 \phi )-\alpha _{\operatorname{ac}} \cos (3 \phi ) \sec (\phi )-3 \alpha _{\operatorname{ac}}\right)}{4 \operatorname{tf1} \left(\mu _2 \cos (2 \phi )-\mu _2-4\right)}-i k \cos (\phi )\right\},\left\{\sigma \to \frac{k^2 \left(-\lambda \alpha _{\operatorname{ac}} \cos ^2(2 \phi )-\alpha _{\operatorname{ac}} \cos (2 \phi )\right)}{\operatorname{tf1} \left(\mu _2 \cos (2 \phi )-\mu _2-4\right)}-\xi \right\}\right\}$