Series coefficients for an infinite sum?
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}}
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.
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}
*)