3
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ What version are you running? I just ran for n = 13 and it took about a minute on 13.3.1 $\endgroup$
    – Edmund
    Dec 18, 2023 at 22:17
  • $\begingroup$ Likewise. Took around 18 minutes for n=17. $\endgroup$ Dec 19, 2023 at 4:39
  • $\begingroup$ 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)? $\endgroup$ Dec 20, 2023 at 9:10

3 Answers 3

5
$\begingroup$

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.

$\endgroup$
5
$\begingroup$

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} *)
$\endgroup$
1
$\begingroup$

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

FunctionProperties`LocallyHolomorphicQ[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.} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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