# Evaluating series expansion is very slow

I need in my work to get series expansion of

(2 E^x x HypergeometricPFQ[{1}, {1/2 + E^-x/4, 1 + E^-x/4}, -(x^2/4)])/Gamma[E^-x/2] + x^(1 - E^-x/2) Sin[x]


up to $$n=20$$. I have tried to run

Normal[Series[(2 E^x x HypergeometricPFQ[{1}, {1/2 + E^-x/4, 1 + E^-x/4}, -(x^2/4)])/Gamma[E^-x/2] + x^(1 - E^-x/2) Sin[x], {x, 0, n}]]


The problem is that time consumed to get the expansion at $$n=13$$ is about 9 minutes, and I waited for hours to evaluate at $$n=20$$ and failed to get a result.

Is there alternative algorithm I may use (or code it myself) to get the result I need?

• What version are you running? I just ran for n = 13 and it took about a minute on 13.3.1 Dec 18, 2023 at 22:17
• Likewise. Took around 18 minutes for n=17. Dec 19, 2023 at 4:39
• Did you write your expressions correctly as E^-x/2 and E^-x/4 or should it be rather E^(-x/2) and E^(-x/4)? Dec 20, 2023 at 9:10

For alternative algorithm and for speed up you can may use: NSeries:

Needs["NumericalCalculus"];

f[x_] := (2 E^x x HypergeometricPFQ[{1}, {1/2 + E^-x/4, 1 + E^-x/4}, -(x^2/4)])/Gamma[E^-x/2] + x^(1 - E^-x/2) Sin[x];

A = NSeries[f[x], {x, 0, 11}, WorkingPrecision -> 100] // Chop // Normal;
N[A /. x -> 1, 30]
(*1.53958275073140462980687669621 - 0.01188149961898482687275300695 I*)

B = Series[f[x], {x, 0, 11}] // Normal;
N[B /. x -> 1, 20]
(*1.5365816625602003563*)

N[f[1], 20]
(*1.5366773917420970181*)


Probably for greater precision you must give more terms.

One can try breaking the expression into pieces:

expr = (2 E^x x HypergeometricPFQ[{1}, {1/2 + E^-x/4, 1 + E^-x/4}, -(x^2/4)])/Gamma[E^-x/2]
+ x^(1 - E^-x/2) Sin[x]
slow[n_] := Normal[Series[expr, {x, 0, n}]]
fast[n_] := Module[{num, den, res},
num = Series[(2 E^x x HypergeometricPFQ[{1}, {1/2 + E^-x/4, 1 + E^-x/4}, -(x^2/4)]), {x, 0, n}];
den = Series[Gamma[E^-x/2], {x, 0, n}];
res = Normal[Series[num/den, {x, 0, n}]] + Normal[Series[x^(1 - E^-x/2) Sin[x], {x, 0, n}]];
res
];


Check:

    Module[{n = 10, f, s},
f = AbsoluteTiming[fast[#]] & /@ Range[n];
s = AbsoluteTiming[slow[#]] & /@ Range[n];
TableForm[Transpose[{Range[n], f[[All, 1]], s[[All, 1]],
Simplify[f[[All, 2]] - s[[All, 2]], Assumptions -> {x > 0}] }],
TableHeadings -> {None, {"n", "Elapsed fast", "Elapsed slow", "Diff"}}
]
]

fast[13]; // AbsoluteTiming
(* {3.74571, Null} *)


The first thing to note, imo, is that the second term of the expression has a branch cut from $$0$$ to $$-\infty$$.

FunctionPropertiesLocallyHolomorphicQ[x^(1 - E^-x/2) Sin[x], x, 0]

(* False *)

FunctionAnalytic[{x^(1 - E^-x/2) Sin[x], Abs[x] < 1/10}, x, Complexes]

(* False *)


It has neither a power nor a Laurent series; rather, its series is a mix of logarithms and square-roots. However, the series is rather quick to compute:

ser2 = Series[x^(1 - E^-x/2) Sin[x], {x, 0, 20}]; // AbsoluteTiming

(* {0.124933, Null} *)

(* x^(11/2) term has a "coefficient" that is a polynomial in Log[x] *)
SeriesCoefficient[ser2, 11/2]

(* (16 + 40 Log[x] + 100 Log[x]^2 - 60 Log[x]^3 + 5 Log[x]^4)/1920 *)


For the first term, one can use NSeries (as @Mariusz did, but with just the first term). Or you can use the Cauchy integral formula:

f = (2 E^x x HypergeometricPFQ[{1}, {1/2 + E^-x/4,
1 + E^-x/4}, -(x^2/4)])/Gamma[E^-x/2];
ser1 = Table[
1/(2 Pi I)
NIntegrate[(f Dt[x, t])/x^(n + 1) /. x -> Exp[I t]/10, {t, 0,
2 Pi}, Method -> "Trapezoidal", WorkingPrecision -> 32,
AccuracyGoal -> 20],
{n, 0, 20}
] // Chop[#, 10^-20] & //
SeriesData[x, 0, #, 0, Length@#, 1] &; // AbsoluteTiming

(* {0.255677, Null} *)


Check the error:

(2 E^x x HypergeometricPFQ[{1}, {1/2 + E^-x/4,
1 + E^-x/4}, -(x^2/4)])/Gamma[E^-x/2] +
x^(1 - E^-x/2) Sin[x] /.
{{x -> 1.32}, {x -> 1/10.32}, {x -> 1/100.32}}

(* {1.5366773917420970181216583991955,
0.14088996522100268533174538168929,
0.012262740086836606574087812170464} *)

Normal[ser1 + ser2] /. x -> {1.32, 1/10.32, 1/100.32}

(* {1.5366774115402366219656862403643,
0.14088996522100268533174556092222,
0.012262740086836606574087812170464} *)

% - %% // N

(* {1.97981*10^-8, 1.79233*10^-25, 0.} *)