I've fiddled with this on and off for a while now, hesitating to decide whether it was worth posting since another answer has already been accepted. The undocumented function, Experimental`OptimizeExpression
, can be used to break down the solutions algebraically into common subexpressions, and it seemed like an approach worth sharing. On the other hand, this equation is essentially equivalent to a quartic polynomial:
Quit[]
eqn = L == (3 W)/2 + (3 Sqrt[4 A^2 Pi^2 + W^2])/2 -
Sqrt[6 A^2 Pi^2 + 3 W^2 + 5 W Sqrt[4 A^2 Pi^2 + W^2]]/Sqrt[2];
sols = W /. Solve[eqn, W];
poly = Collect[5 Times @@ (W - sols) // Expand // Simplify, W]
(*
L^4 - 24 A^2 L^2 π^2 + 36 A^4 π^4 + (-6 L^3 + 12 A^2 L π^2) W +
(6 L^2 + 20 A^2 π^2) W^2 - 6 L W^3 + 5 W^4
*)
The solutions turn out to resemble the standard quartic formula. Indeed, it might be just as easy to build up the solutions from the formula by hand from the polynomial derived below (perhaps using Mathematica to keep track of the algebraic steps). The expression will still be rather complicated due to the nature of the quartic formula and the coefficients, even if they appear to have some symmetry.
Different purposes might be fulfilled by being able to write expressions for the roots. Presumably the general goal would be to illuminate elements of the problem, but all we are given is an algebraic equation. In some cases an analytic approach might be more appropriate than an algebraic approach, or vice versa.
###Auxiliary functions###
Set up the equation and solutions; the solutions may be used to construct a polynomial with the same roots (the coefficient 5
is discovered by inspection).
The function Experimental`OptimizeExpression
returns an expression of the form
Experimental`OptimizedExpression[Block[{vars}, var1 = val1; <>; varn = valn; expr]]
The variables have the form Compile`$nnn
, where nnn
represents a serial number, and represent subexpressions; if a subexpression appears more than once, it will be represented by the same variable. Unfortunately, the serial number increases throughout a session and the starting number depends on what evaluations have been done. (Execute Quit[]
above to get the same numbering -- unless the numbering is version/system dependent.) The expression expr
represents the optimized expression in terms of the variables. Preceding it, the variables are initialized. From this we can construct auxiliary functions for exploring the expression.
Experimental`OptimizeExpression
takes an OptimizationLevel
option and may be set to 0
, 1
, or 2
. Subexpressions are stored in the returned Experimental`OptimizeExpression
in a held Block
in the form Compile`$nnn = expr
. We can turn these Set
expressions into Rule
expressions that can be used to expand a given subexpression in terms of the next level of subexpressions. One can also replace the last expression inside Block
with an arbitrary expression in terms of Compile`$nnn
variable and it will be evaluated in terms of the subexpressions of the optimized expression.
optexpr = Experimental`OptimizeExpression[sols, OptimizationLevel -> 2];
(* get the initialization in terms of Rule *)
optrules = Most[(optexpr /. {Set -> Rule, CompoundExpression -> List})[[1, 2]]];
(* convert between Compile`$nnn symbol and the number nnn *)
compileSym = ToExpression["Compile`$" <> ToString[#]] &;
compileSymNo = ToExpression@StringDrop[SymbolName[#], 1] &;
(* get the range of the serial numbers of the variables *)
{minCompileNumber,
maxCompileNumber} =
Through[{Min, Max}[
Cases[optexpr,
x_Symbol /; Context[x] === "Compile`" :> compileSymNo[x](*ToExpression@
StringDrop[SymbolName[x],1]*), Infinity]]];
(* evaluate an expression in terms of Compile`$nnn variables *)
evalCompileExpr = ReplacePart[Function @@ optexpr, {1, -1, -1} :> Slot[1]];
The following shows the structure of the four roots determined by Experimental`OptimizeExpression
. It reveals the typical structure of the roots of a quartic equation.
optoutput = optexpr[[1, -1, -1]]
(*
{Compile`$1 + Compile`$57 - Compile`$71/2,
Compile`$1 + Compile`$57 + Compile`$71/2,
Compile`$1 + Compile`$76 - Compile`$79/2,
Compile`$1 + Compile`$76 + Compile`$79/2}
*)
###Using the functions###
Executing expr /. optrules
will expand the variables one step. What steps to take requires some judgment. If we expand optoutput
, we see the structure of the quartic formula begin to unfold.
optoutput /. optrules
(*
{(3 L)/10 - Compile`$56/2 - Sqrt[Compile`$70]/2, (3 L)/10 - Compile`$56/2 +
Sqrt[Compile`$70]/2, (3 L)/10 + Compile`$56/2 - Sqrt[Compile`$78]/2, (3 L)/
10 + Compile`$56/2 + Sqrt[Compile`$78]/2}
*)
{Compile`$70, Compile`$78} /. optrules
(*
{Compile`$58 + Compile`$59 + Compile`$60 + Compile`$61 + Compile`$69,
Compile`$58 + Compile`$59 + Compile`$60 + Compile`$61 + Compile`$77}
*)
{Compile`$69, Compile`$77} /. optrules
(*
{-((Compile`$67 Compile`$68)/4), (Compile`$67 Compile`$68)/4}
*)
The above calculations show we can write the solutions in the form
W -> 3 L / 10 ± d1 ± Sqrt[d2 ± d3] / 2
where
d1 = Compile`$56 / 2
d2 = Compile`$58 + Compile`$59 + Compile`$60 + Compile`$61
d3 = (Compile`$67 Compile`$68) / 4
and the first and third ± signs agree.
One can examine the actual expressions with evalCompileExpr
, but, as I said, without the context of the problem, it's hard to see anything important lurking in the expressions.
We can break down the solution into bit-size pieces -- well, whiteboard-size pieces. One can see there is a large cube root that is repeated and larger square root. We can get them and replace with them with new variables Q
and R
as follows:
rootTerms = {Compile`$56/2,
Compile`$58 + Compile`$59 + Compile`$60 + Compile`$61, (Compile`$67 Compile`$68)/4};
simpRT = Simplify /@ evalCompileExpr[rootTerms];
q = First@Cases[simpRT, Power[Except[_?NumberQ], 1/3], Infinity];
Q -> q
cbrt = {q -> Q, 1/q -> 1/Q};
r = First@Cases[simpRT /. cbrt, Power[Except[_?NumberQ], 1/2], Infinity];
R -> r
sqrt = {r -> R, 1/r -> 1/R};
Thread[{d1, d2, d3} -> simpRT /. cbrt /. sqrt]
(*
Q -> (135 L^6 + 225 A^2 L^4 π^2 + 11052 A^4 L^2 π^4 - 15200 A^6 π^6 +
3 Sqrt[3] Sqrt[(13 L^2 + 24 A^2 π^2)^2 (4 L^8 + 767 A^4 L^4 π^4 - 2000 A^8 π^8)])^(1/3)
R -> Sqrt[27 L^2 - 20 (3 L^2 + 10 A^2 π^2) -
(10 (3 L^4 + 246 A^2 L^2 π^2 - 640 A^4 π^4))/Q + 10 Q]
{d1 -> R/(10 Sqrt[3]),
d2 -> 1/75 (54 L^2 - 40 (3 L^2 + 10 A^2 π^2) +
(10 (3 L^4 + 246 A^2 L^2 π^2 - 640 A^4 π^4))/Q - 10 Q),
d3 -> (6 Sqrt[3] (29 L^3 - 200 A^2 L π^2))/(25 R)}
*)
These together with the expression of W
above in terms of d1
, d2
, and d3
present the complete solution. Aside from some minor simplifications, one can see that d2
may be written
d2 -> 1/75 (81 L^2 - 60 (3 L^2 + 10 A^2 π^2) - R^2) // Factor
(*
d2 -> 1/75 (-99 L^2 - 600 A^2 π^2 - R^2)
*)
###Homogenization###
Here is another way to look at the solutions. It is not amenable to written presentation but it is a nice way to look at the problem. With a change of variables, we can transform the equation into a homogeneous polynomial poly0
of three variables. Dilations then act on the solution set and one dimension can be factored out. In other words all solutions may be obtained from scaling a given cross-section of the surface poly0 == 0
.
Here are transformations for converting between poly
and poly0
.
homogenize = {A :> Sqrt[α]/(Sqrt[2] Pi), L -> λ/Sqrt[α], W -> Ω/Sqrt[α]};
dehomogenize = First@Solve[{A, L, W} == ({A, L, W} /. homogenize), {α, λ, Ω}]
(*
{α -> 2 A^2 π^2, λ -> Sqrt[2] A L π, Ω -> Sqrt[2] A π W}
*)
eqn0 = eqn /. homogenize // Simplify;
sols0 = Ω /. Solve[eqn0, Ω];
poly0 = Collect[5 Times @@ (Ω - sols0) // Expand // Simplify, Ω]
poly == ((W/Ω)^4 poly0 /. dehomogenize) // Expand
(*
9 α^4 - 12 α^2 λ^2 + λ^4 + (6 α^2 λ - 6 λ^3) Ω + (10 α^2 + 6 λ^2) Ω^2 - 6 λ Ω^3 + 5 Ω^4
True
*)
Any solution (for Ω
or A π W
) to poly0 == 0
may be obtained by dilation (scaling) of a plane section of the surface poly0 == 0
. Below the relationships of each of two sections to the solution set are shown, with the mesh lines showing the dilation of the boundary curve.
sectλ = Show[
ContourPlot3D[
poly0 == 0, {α, -1.2, 1.2}, {λ, -1, 1}, {Ω, -1, 1},
MeshFunctions -> {ArcTan[#1, #2] &, #2 &},
PlotPoints -> 20, AxesLabel -> Automatic],
ParametricPlot3D[
Thread[{α, λ, sols0}] /. λ -> 1 // Evaluate, {α, -1.2, 1.2},
PlotPoints -> 100,
PlotStyle -> (Directive[Thickness[0.01], #] & /@ {Red, Blue, Magenta, Darker@Green})]
];
sectα = Show[
ContourPlot3D[
poly0 == 0, {α, -1, 1}, {λ, -4, 4}, {Ω, -3, 3},
MeshFunctions -> {ArcTan[#1, #2] &, #1 &}, PlotPoints -> 20,
AxesLabel -> Automatic],
ParametricPlot3D[
Thread[{α, λ, sols0}] /. α -> 1 // Evaluate, {λ, -4, 4},
PlotPoints -> 100,
PlotStyle -> (Directive[Thickness[0.01], #] & /@ {Red, Blue, Magenta, Darker@Green})],
BoxRatios -> {1, 4, 3}
];
GraphicsRow[{sectλ, sectα}]
Solve
). $\endgroup$Simplify
orFullSimplify
? Some assumptions might help with those, too. $\endgroup$