3
$\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.

$\endgroup$
  • 1
    $\begingroup$ Series[ToRadicals@Root[j^2 - 256 #1^2 + 1024 #1^3 &, 2], {j, 0, 4}]? $\endgroup$ – kglr Mar 2 '15 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 '15 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 '15 at 9:34
  • $\begingroup$ Thank you @Alexei; just posted the comment as an answer. $\endgroup$ – kglr Mar 4 '15 at 9:40
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!

$\endgroup$
3
$\begingroup$

Perhaps

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

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Yes, thank you! $\endgroup$ – Alexei Boulbitch Mar 4 '15 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 '17 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.

$\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.

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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