Series coefficients for an infinite sum?

Question

I'm working through Carl Bender's Mathematical Physics lectures on YouTube (which are great fun), and I'd like Mathematica's help solving terms in the perturbation series. It would be convenient if expressions like

SeriesCoefficient[Sum[a[n] b^n, {n, 0, ∞}], {b, 0, 3}]

gave a[3] as a result. Replacing $\infty$ with an integer greater than 3 is fine, but easy to lose track of in my notes. Is it possible to add an assumption that the series converges, and would that solve my problem? Other solutions?

Edit

For an example of how I hope to use this functionality, consider solving $x^5 + x = 1$ by inserting the perturbation parameter $b$ in front of $x$ (but not $x^5$). I want to apply SeriesCoefficient to the result of

x^5 + b x - 1 /. x -> Sum[a[n] b^n, {n, 0, ∞}]

and build up a list of coefficients that equal zero.

Answer 1

How funny - I was just doing this the other day (different equation, with a series in x^(1/2)). Here is what I came up with (actually, I'm following Isaac Newton's method, I believe):

Clear[sc];
(* series coefficient of x in f[b, x] == 0 *)
sc[f_, n_] := Nest[
  Append[#, a /. First @ Solve[
       SeriesCoefficient[
         f[b, #.b^(Range @ Length @ # - 1) + a b^(Length @ #)],
         {b, 0, Length @ #}] == 0,
       a]] &,
  {a /. First @ Solve[SeriesCoefficient[f[b, a], {b, 0, 0}] == 0, a]},
  n]

OP's example:

ClearAll[f];
f[b_, x_] := x^5 + b x - 1

sc[f, 5]
(*
   {1, -(1/5), -(1/25), -(1/125), 0, 21/15625}
*)

sc[f, 25] // Total // N
(*
   0.754878
*)

NSolve[f[1, x] == 0, Reals]
(*
   {{x -> 0.754878}}
*)

---

There are imaginary roots too:

NRoots[f[1, x] == 0, x]
(*
  x == -0.877439 - 0.744862 I || x == -0.877439 + 0.744862 I ||
  x == 0.5 - 0.866025 I || x == 0.5 + 0.866025 I || x == 0.754878
*)

To get all the roots in OP's example:

Clear[sc];

sc[f_, n_] := Nest[
  Map[Join[#,
    {a} /. First@Solve[SeriesCoefficient[f[b, #.b^(Range@Length@# - 1) + a b^(Length@#)],
    {b, 0, Length@#}] == 0, a]] &, #] &,
  Transpose[{a /. Solve[SeriesCoefficient[f[b, a], {b, 0, 0}] == 0, a]}],
  n]

Total[sc[f,5],{2}]//N
(*
  {0.753344, -0.877796 - 0.745447 I, 0.501124 + 0.865898 I,
  0.501124 - 0.865898 I, -0.877796 + 0.745447 I}
*)

Answer 2

The new in M12 function AsymptoticSolve can be used to find the perturbation expansions:

AsymptoticSolve[x^5 + b x == 1, {x, 1}, {b, 0, 5}]

{{x -> 1 - b/5 - b^2/25 - b^3/125 + (21 b^5)/15625}}

Or, perturbations around all 5 roots can be found with:

AsymptoticSolve[x^5 + b x == 1, x, {b, 0, 5}]

{{x -> 1 - b/5 - b^2/25 - b^3/125 + (21 b^5)/15625}, {x -> -(-1)^(1/5) - 1/5 (-1)^(2/5) b + 1/25 (-1)^(3/5) b^2 - 1/125 (-1)^(4/5) b^3 - ( 21 (-1)^(1/5) b^5)/15625}, {x -> (-1)^(2/5) - 1/5 (-1)^(4/5) b + 1/25 (-1)^(1/5) b^2 + 1/125 (-1)^(3/5) b^3 + (21 (-1)^(2/5) b^5)/ 15625}, {x -> -(-1)^(3/5) + 1/5 (-1)^(1/5) b - 1/25 (-1)^(4/5) b^2 - 1/125 (-1)^(2/5) b^3 - (21 (-1)^(3/5) b^5)/15625}, {x -> (-1)^(4/5) + 1/5 (-1)^(3/5) b - 1/25 (-1)^(2/5) b^2 + 1/125 (-1)^(1/5) b^3 + ( 21 (-1)^(4/5) b^5)/15625}}

Answer 3

You could always build a pattern matcher that will recover the coefficients you want:

seriescoefficient[
  Sum[coeff_ b_^n_, {n_, min_, max_}], {b_, b0_, m_}] := coeff /. {n -> m}

will evaluate

seriescoefficient[Sum[a[n] b^n, {n, 0, ∞}], {b, 0, 3}]

to a[3]. This may or may not work, or require different amounts of care for conditions, depending on how narrow a class of inputs you want to handle. The code above, for example, will return the same input even if the Sum's iterator starts from n==5; this can of course be fixed but the question is where do you stop fixing potential trouble inputs.

You can also add this behaviour to the old SeriesCoefficient function by running

Unprotect[SeriesCoefficient]

and then the definition above (be sure to Protect it when you're done, though).

---

The approach above is useful whenever you can express whatever you want to simplify as a single series. However, Mathematica will not willingly convert powers and products of series into more complicated single series. Another approach that can work is to simply take the appropriate derivative.

The problem with the infinite sums is that they are not always great at substituting in b->0, particularly when the power b^0 is present. The function

substitutor[series_] := (series /. Sum -> sum) /. {
    sum[b_^(exp_ + n_) coeff_, {n_, 0, ∞}] -> replaceall[coeff, n -> -exp],
    sum[b_^n_ coeff_, {n_, 0, ∞}] -> replaceall[coeff, n -> 0]
    } /. replaceall -> ReplaceAll

is pretty messy, but it does the job of taking a sum of the form Sum[b_^(-n0+n) patt_, {n, 0, \[Infinity]}] and returning its value when b==0. With this, for example,

substitutor[D[
  ((x^5 + x - 1) /. x -> Sum[a[n] b^n, {n, 0, ∞}])
  , {b, 3}]]

returns the same as what

3! SeriesCoefficient[((x^5 + x - 1) /. x -> Sum[a[n] b^n, {n, 0, 15}]), {b, 0, 3}]

does. Again, this may or may not be what you need, but it can work.

Answer 4

The (accepted) answer given by episanty gave me another idea I want to share. While SeriesCoefficient leaves the infinite sum untouched, Derivative does not, so we can get the coefficients manually. The issue episanty points out is that the lower limit of the sum errors out at $0^0$. We can prevent that behavior by modifying Sum (eek!) and everything works out nicely. I don't know what other problems it will create, but here is one sufficient downvalue that is also pretty narrowly defined:

Sum[0^k_ a_, {n_, 0, DirectedInfinity[1]}] := a /. First@Solve[k == 0, n]

You have to Unprotect[Sum] first, of course, and then protect it again.

Now, to get equations for the first 4 coefficients:

f[b_] := x^5 + b x - 1 /. x -> Sum[a[n] b^n, {n, 0, ∞}]

Simplify@Table[1/k! Derivative[k][f][0] == 0, {k, 0, 3}]

To restore the default Sum, execute DownValues[Sum] = {}.