4
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Series[ToRadicals@Root[j^2 - 256 #1^2 + 1024 #1^3 &, 2], {j, 0, 4}]? $\endgroup$
    – kglr
    Mar 2, 2015 at 21:53
  • $\begingroup$ 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][[2]]], {J, 0, 4}] $\endgroup$
    – Dr. belisarius
    Mar 2, 2015 at 22:03
  • $\begingroup$ @kguler I like you solution, but why do not you formulate your comment as a regular answer? $\endgroup$
    – Alexei Boulbitch
    Mar 4, 2015 at 9:34
  • $\begingroup$ Thank you @Alexei; just posted the comment as an answer. $\endgroup$
    – kglr
    Mar 4, 2015 at 9:40

5 Answers 5

Reset to default
3
$\begingroup$

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}], #[[2]] <= 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:

enter image description here

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: enter image description here

Have fun!

Improve this answer
$\endgroup$
3
$\begingroup$

Perhaps

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

enter image description here

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Yes, thank you! $\endgroup$
    – Alexei Boulbitch
    Mar 4, 2015 at 10:02
  • 2
    $\begingroup$ Do you know any way of finding the series in case the root has no radical representation? $\endgroup$
    – Thomas Ahle
    May 20, 2017 at 0:50
2
$\begingroup$

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.

Improve this answer
$\endgroup$
1
$\begingroup$

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:

Visualization of the three branches.

Improve this answer
$\endgroup$
1
$\begingroup$

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}}

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.