# Series expression of a Root object

I was wondering if it is possible to get Mathematica to return a series approximation of a Root object. Example: I want a series representation of x in terms of J for the equation

16 x^3 - 4 x^2 + J^2/64 == 0


When I use Solve I get

Root[J^2 - 256 #1^2 + 1024 #1^3 &, 2]


for some domain of J. So I try

Series[Root[J^2 - 256 #1^2 + 1024 #1^3 &, 2], {J, 0, 4}]


which returns

Series::nmer: "Root[J^2 - 256 #1^2 + 1024 #1^3 &, 2] is not a meromorphic function of J at 0."

which makes sense because the object is only defined for a certain domain of J. I was wondering if there is a way to get around this and force Mathematica to output a series representation for that root as a function of J in that domain.

• Series[ToRadicals@Root[j^2 - 256 #1^2 + 1024 #1^3 &, 2], {j, 0, 4}]?
– kglr
Mar 2, 2015 at 21:53
• Or expanding the above comment by kguler :Series[ToRadicals@ First[x /. First@Solve[16 x^3 - 4 x^2 + J^2/64 == 0, x, Reals][]], {J, 0, 4}] Mar 2, 2015 at 22:03
• @kguler I like you solution, but why do not you formulate your comment as a regular answer? Mar 4, 2015 at 9:34
• Thank you @Alexei; just posted the comment as an answer.
– kglr
Mar 4, 2015 at 9:40

There are excellent answers of kguler and belisarius. Not to compete with them, but to offer a different view, you might operate within a different paradigm. Here it is:

Step 1: Let us first find the solution of your equation for J:

 sl = Solve[16 x^3 - 4 x^2 + J^2/64 == 0, J]

(*   {{J -> -16 Sqrt[x^2 - 4 x^3]}, {J -> 16 Sqrt[x^2 - 4 x^3]}}  *)


There are two, and let us operate on the positive one. The negative may be treated analogously.

Step 2: Let us create a table entitled "lst" out of this second solution with the structure {J,x}. It is easily seen that x=x(J)is a two-valued function and I only select the smaller branch:

lst = Select[
Table[{sl[[2, 1, 2]], x}, {x, 0, 0.25, 0.01}], #[] <= 0.17 &]

(*  {{0., 0.}, {0.156767, 0.01}, {0.306933, 0.02}, {0.45028,
0.03}, {0.58657, 0.04}, {0.715542, 0.05}, {0.836909,
0.06}, {0.950352, 0.07}, {1.05552, 0.08}, {1.152, 0.09}, {1.23935,
0.1}, {1.31706, 0.11}, {1.38453, 0.12}, {1.44107, 0.13}, {1.48585,
0.14}, {1.51789, 0.15}, {1.536, 0.16}, {1.53866, 0.17}}  *)


The upper one can be selected analogously.

Step 3: Let us now approximate (that is, fit) the obtained table with a polynomial:

 model1 = a + b*J + c*J^2 + d*J^3;
ff1 = FindFit[lst, model1, {a, b, c, d}, J]
model2 = a + b*J + c*J^2 + d*J^3 + e*J^4;
ff2 = FindFit[lst, model2, {a, b, c, d, e}, J]

(*  {a -> -0.00294399, b -> 0.106524, c -> -0.0890299, d -> 0.0569192}

{a -> 0.00137571, b -> 0.0172211, c -> 0.186767, d -> -0.219333,
e -> 0.0872711}   *)


The first solution is for the first model (that is, the expansion up to the cube), while the second is for the second one (the expansion up to J^4). Let us draw them to see the difference:

 Show[{
ListPlot[lst, AxesLabel ->  {"J", "x"}],
Plot[{model1 /. ff1, model2 /. ff2}, {J, 0, 1.5},
PlotStyle -> {Red, Green}]
}]


It returns the following plot: The red here shows the fitting with the first model, while green - with the second one. Since this fitting yields a polynomial this is though not equal to, but not really too far away from what kguler and belisarius proposed. The real difference might come, if the expansion is not a final aim, but a step in some more general problem. Then one might think about a simple and exact enough approximation. Just to give a simple example, have a look at the following two models with a smaller amount of parameters:

 model3 = (a*J)/(1 + d*J^2);
ff3 = FindFit[lst, model3, {a, d}, J]
model4 = (a*J + b*J^2)/(1 + d*J^2);
ff4 = FindFit[lst, model4, {a, b, d}, J]

(*  {a -> 0.0596691, d -> -0.176778}
{a -> 0.0807892, b -> -0.0277843, d -> -0.269206}  *)


The visualization is below: Have fun!

Perhaps

Series[ToRadicals@Root[j^2 - 256 #1^2 + 1024 #1^3 &, 2], {j, 0, 4}] • Yes, thank you! Mar 4, 2015 at 10:02
• Do you know any way of finding the series in case the root has no radical representation? May 20, 2017 at 0:50

In version 11, ToRadicals[] is no longer necessary:

Series[Root[J^2 - 256 #1^2 + 1024 #1^3 &, 2], {J, 0, 5}]
J/16 + J^2/128 + (5 J^3)/2048 + J^4/1024 + (231 J^5)/524288 + O[J]^6


altho the warning Root::sbr is also given.

Alternatively, you may reverse the order and expand in J before solving the equation:

eq = 16 x^3 - 4 x^2 + J^2/64;
sol = CoefficientList[Series[eq /. x -> x[J], {J, 0, 5}], J] == 0 // LogicalExpand // Solve;
Normal@Series[x[J], {J, 0, 4}] /. sol


Gives you approximations for three branches:

{-(J/16) + J^2/128 - (5 J^3)/2048 + J^4/1024, J/16 + J^2/128 + (5 J^3)/2048 + J^4/1024, 1/4 - J^2/64 - J^4/512}


This is what these approximations look like to 4th order: The new in M12 function AsymptoticSolve can be used to find the series:

AsymptoticSolve[16 x^3 - 4 x^2 + J^2/64 == 0, x, {J, 0, 4}]


{{x -> -(J/16) + J^2/128 - (5 J^3)/2048 + J^4/1024}, {x -> J/16 + J^2/128 + (5 J^3)/2048 + J^4/1024}, {x -> 1/4 - J^2/64 - J^4/512}}