# Solve of cubic and quartic equations too slow

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

• I find that sol = Solve[eq == 0, σ] takes just seconds, although the result is rather large. By the way, avoid using Subscript. – bbgodfrey Dec 9 '15 at 0:20
• Furthermore, solving the fourth order equation takes only a minute or so. How much memory do you have on your computer? – bbgodfrey Dec 9 '15 at 0:32
• sol = Solve[eq == 0, σ] took 2 minutes on my computer. The method I presented below takes around 3 seconds. – JungHwan Min Dec 9 '15 at 1:01
• bbgodfrey, what PC do you have? I have laptop with 8GB RAM, quadcore intel i7. Using device manager I found out, that during the calculations ~20% of CPU and ~40% of RAM is used, so I doubt, that the issue is in memory. Thanks for the edit, btw – Mikhail Genkin Dec 9 '15 at 1:26
• @MikhailGenkin I have an essentially identical computer running 10.3.0 for Microsoft Windows (64-bit) (October 9, 2015). The calculations complete so quickly that I cannot accurately determine CPU and memory usage with Task Manager for the cubic case. The quartic case takes less than a minute with 17% CPU usage and about 200 MB memory usage. – bbgodfrey Dec 9 '15 at 1:36

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.

• Thanks a lot, you method really works fast on my laptop! Can you briefly explain, why? But this thing definitely works for me. I will accept your answer it 2 days. – Mikhail Genkin Dec 9 '15 at 1:23
• Nevermind, I got it. The reason I didn't figure it out by myself, is because I thought, that Solve[eq==0,s] works like this. It seems that @bbgodfrey has a new version, which is doing it by default. (He wrote, that his mathematica is solving it very fast). I assume, that you too have an older version of Mathematica, like I do. My is 10.0.1.0 – Mikhail Genkin Dec 9 '15 at 1:43
• Nice but might be overly complicated. To see what I mean, compare constant terms with what you get from Solve[(eq /. k -> 0) == 0, \[Sigma]]. Much smaller and in agreement with the larger result after some random substitution. – Daniel Lichtblau Dec 9 '15 at 15:12

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];
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.

• Thanks for the interesting method. Your code seems to be incorrect: my result is different. You never use variable derivs – Mikhail Genkin Dec 9 '15 at 19:04
• (1) Thanks for spotting that mistake. I had some earlier definitions and decided to change notation; what got posted was a bit of a mixture. I believe I have corrected it now. – Daniel Lichtblau Dec 9 '15 at 19:41
• (2) It is possible that the solution sets are equivalent. One way to check is to substitute values in for most or all parameters (other than k of course) and see if the sets of three numeric series agree. – Daniel Lichtblau Dec 9 '15 at 19:42
• I checked it against direct method. Solution are equivalent, as should be. I also got the results really fast for 4th order equation, which is not the case for direct method. You helped me a lot. – Mikhail Genkin Dec 9 '15 at 20:10
• I still accepted other answer. Even though your answer is usefull, the other one directly answers my question, whether yours is more workaround. Thanks again – Mikhail Genkin Dec 11 '15 at 17:39

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[SystemAsymptoticSolve[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\}$$