# Series coefficients for an infinite sum?

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.

• I think that SeriesCoefficient will only work with closed expressions. If the sum inside fails to evaluate, so will the SeriesCoefficient. But I may be wrong. Feb 27, 2014 at 13:53
• @m_goldberg How did you add the infinity symbol in the code block? And why would so magical a being not fix the ugly -> too?
– Ian
Feb 27, 2014 at 16:51
• Additional useful buttons for our M.SE editor. The reason the -> doesn't change is presumably that there isn't an ASCII character for it that MM will interpret correctly. Feb 27, 2014 at 17:01
• In this case I used a browser plug-in written by halirutan, but for inserting special characters, you just paste them in from a unicode character palette. Feb 27, 2014 at 18:53

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


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

• Somehow, I am not able to check this... :/ Apr 14, 2019 at 1:34
• @MichaelE2 At least, I think you can use AsymptoticSolve on the cloud for now. Apr 14, 2019 at 2:31
• Thanks! That works. Apr 14, 2019 at 3:13

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.

• Do you think this could be expanded to address arbitrary expressions involving the infinite sum? How I really want to use it ... well, let me go back and put that in the question.
– Ian
Feb 27, 2014 at 13:53

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] = {}.