How to make faster of a code containing more iteration points?

Question

I am trying to evaluate a code of the Power Series Solution method to get the results. Also, the time of evaluation depends on the parameters I am using. When the iteration point is 50, 60, or 80, it takes less time(10 mins to 1 hour) than when the iteration point is 100 or 120(like two to three days). So is there any code that makes it faster to get the result? Here is an example. When I am giving the iteration value 100, it takes two days whereas in the case of 40, 50, or 80 it takes 15-20 mins, and for smaller values of $qE$ it takes very large time.

g[x_] = (1 + 
      1/(4 a^4) (8 a^2 (a z - Log[1 + a z]) + (-2 a^2 + qE^2 + 
            qM^2) (a z (6 - a z) + 2 (-3 - 2 a z + a^2 z^2) Log[1 + a z] + 
            2 Log[1 + a z]^2) + 
         2 a^2 (a z (4 - a z) + 2 Log[z] (a z (-2 + a z) + 2 Log[1 + a z]) + 
            4 PolyLog[2, -a z])) - 
      1/(a zh (-2 + a zh) + 
        2 Log[1 + a zh]) (a z (-2 + a z) + 2 Log[1 + a z]) (1 + 
         1/(4 a^4) (8 a^2 (a zh - Log[1 + a zh]) + (-2 a^2 + qE^2 + 
               qM^2) (a zh (6 - a zh) + 
               2 (-3 - 2 a zh + a^2 zh^2) Log[1 + a zh] + 
               2 Log[1 + a zh]^2) + 
            2 a^2 (a zh (4 - a zh) + 
               2 Log[zh] (a zh (-2 + a zh) + 2 Log[1 + a zh]) + 
               4 PolyLog[2, -a zh])))) /. z -> x /. zh -> xh;

A[x_] = -Log[1 + a x];
F[x_] = g[x]/x^2;
G[x_] = g[x]/(x^2 Exp[2 A[x]])

(* (1/(x^2))(1 + a x)^2 (1 + (
   8 a^2 (a x - Log[1 + a x]) + (-2 a^2 + qE^2 + qM^2) (a x (6 - a x) + 
       2 (-3 - 2 a x + a^2 x^2) Log[1 + a x] + 2 Log[1 + a x]^2) + 
    2 a^2 (a x (4 - a x) + 2 Log[x] (a x (-2 + a x) + 2 Log[1 + a x]) + 
       4 PolyLog[2, -a x]))/(
   4 a^4) - ((a x (-2 + a x) + 2 Log[1 + a x]) (1 + (
      8 a^2 (a xh - Log[1 + a xh]) + (-2 a^2 + qE^2 + 
          qM^2) (a xh (6 - a xh) + 2 (-3 - 2 a xh + a^2 xh^2) Log[1 + a xh] + 
          2 Log[1 + a xh]^2) + 
       2 a^2 (a xh (4 - a xh) + 
          2 Log[xh] (a xh (-2 + a xh) + 2 Log[1 + a xh]) + 
          4 PolyLog[2, -a xh]))/(4 a^4)))/(a xh (-2 + a xh) + 2 Log[1 + a xh])
   ) *)

Scoeff = -((x^4 G[x])/(x - xh));
Scoeffa0 = Series[Scoeff, {a, 0, 0}] // Normal

(* -((x^2 (-4 x^2 - 4 x xh - 4 xh^2 - 4 x^2 xh^2 + qE^2 x^3 xh^3 + 
    qM^2 x^3 xh^3))/(4 xh^3)) *)

Clear[xh]
Scoeffa0 /. qE -> 0 /. qM -> 0 // Expand

(* x^4/xh^3 + x^3/xh^2 + x^2/xh + x^4/xh *)

tcoeff = -2 x^3 G[x] - 2 I x^2 ω Sqrt[G[x]/F[x]] - (
   x^4 G[x] Derivative[1][F][x])/(2 F[x]) - 1/2 x^4 Derivative[1][G][x];

tcoeffa0 = Series[tcoeff, {a, 0, 0}] // Normal // Expand

(* -2 x^3 - qE^2 x^5 - qM^2 x^5 + (3 x^4)/xh^3 + (3 x^4)/xh + 3/4 qE^2 x^4 xh + 3/4 qM^2 x^4 xh - 2 I x^2 ω *)

ucoeff = L x^3 + L^2 x^3 - L x^2 xh - L^2 x^2 xh - (
    x^4 G[x] Derivative[1][F][x])/(2 F[x]) + (
    x^3 xh G[x] Derivative[1][F][x])/(2 F[x]) - 1/2 x^4 Derivative[1][G][x] + 
    1/2 x^3 xh Derivative[1][G][x] // Simplify;

ucoeffa0 = Series[ucoeff, {a, 0, 0}] // Normal // Expand;

QNM Calculation

Clear[qE, qM, s, t, u, Sn, Tn, Un, P, aa, ex, expressao, rh, xh]

L = 0;
qM = 0;
qE = 5;
(*a=0.01;*)
rh = 2; xh = 1/rh;
s[x_] = Scoeffa0 // Simplify
t[x_, ω_] = tcoeffa0 // Simplify
u[x_] = ucoeffa0 // Simplify

(* 
1/4 x^2 (8 + 16 x + 40 x^2 - 25 x^3)
-2 x^3 + (315 x^4)/8 - 25 x^5 - 2 I x^2 ω
-1 + 2 x - (105 x^3)/16 + (155 x^4)/8 - (25 x^5)/2 *)

Sn = Function[n, SeriesCoefficient[Series[s[x], {x, xh, n}], n]];
Tn = Function[n, SeriesCoefficient[Series[t[x, ω], {x, xh, n}], n]];
Un = Function[n, SeriesCoefficient[Series[u[x], {x, xh, n}], n]];

P = Function[{n, ω}, n (n - 1) Sn[0] + n Tn[0]];

Clear[aa, Nit, ex]
aa[0] = 1;
Nit = 100;
 aa[n_] := 
  aa[n] = Simplify[-1/
      P[n, ω] Sum[(k (k - 1) Sn[n - k] + k Tn[n - k] + Un[n - k]) aa[
        k], {k, 0, n - 1}]];

ex[n_] := ex[n] = Sum[aa[i]*(-xh)^i, {i, 0, n}];
expressao := Sum[aa[i]*(-xh)^(i), {i, 0, Nit}]

QNMs = FindRoot[expressao == 0, {ω, 100 - 200*I}, AccuracyGoal -> 10][[1]][[2]]
wroots = SolveValues[expressao == 0, ω] // N```

Answer 1

The code runs for me in 3-4 mins with Nit=100. The only thing I changed was the WorkingPrecision in FindRoot and SolveValues (since you want numerical approximations anyways). I also removed ex because it's unused.

(*all the same code as yours*)

(*...until here*)
ClearAll[expressao]
AbsoluteTiming[
 expressao := Sum[aa[i]*(-xh)^(i), {i, 0, Nit}];
 QNMs = FindRoot[expressao == 0, {\[Omega], 100 - 200*I}, 
     WorkingPrecision -> 64][[1]][[2]];
 wroots = 
  SolveValues[expressao == 0, \[Omega], WorkingPrecision -> 128];
 ]

{221.373, Null}

And verifying the roots:

Chop[expressao /. \[Omega] -> QNMs]
Chop[expressao /. \[Omega] -> wroots]

0

{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Answer 2

This might not fully solve your problem, but I believe it's worth the try:

All the functions you are using here are holomorphic functions (in fact, I think I only saw hyper-geometric ones). This means that you can re-write everything in the form: $$S(x)=\sum_n a_n x^n$$ and in fact, even in the form $$S(x)=\sum_{n_1,n_2,n_3} a_{n_1,n_2,n_3} x^{n_1} a^{n_2} h^{n_3}$$ if it turns out you want to try several values of $a$ and $h$.

The nice thing about (actually) computing everything in this form is that the coefficients $a_n$ are always the same between each execution. This gives you the luxury of storing them somewhere, so that you only have to compute them once (after that, the complexity is only that of reading data in a file, (i.e. $\mathcal{O}(1)$). As for the changing part of the computation (i.e. the "$x^n$") it can be recursively stored at each step so that you are computing it also in $\mathcal{O}(1)$, bringing your total computation of a quick $\mathcal{O}(N)$ (where $N$ is the number of terms in your series). Worst case scenario (if you go to higher $N$s each time, the fact that your series is hypergeometric means each terms can be computed using the previous, meaning the total computation is still $\mathcal{O}(N)$ (provided you are computing in floats, if you are computing in exact values, your numerator and denominator of $a_n$ are going to blow up like as products of factorials, and so the reduction of the fraction is going to tremendously slow your computer down). But in any case, you will only have to do that once for each term since you end up storing it.

The bad thing about this otherwise ideal method (in spite of all the pain you are going to go through while figuring out what $a_n$ is) is that $S$ has a radius of convergence (or two, if you are doing a Laurent expansion). This means that if you are computing for $x$ sometimes inside the convergence area and sometimes not then you need to use another series. And there's no telling how many of these you'll need (potentially: an infinite number; worst case: one for each of the values you want to try) and that quite defeats the purpose of using series in the first place.

To conclude with the way to proceed:

• first consider the full range over which you want to evaluate.
• try to see if you can find singularities in your function; they are a good indicator of the size of your convergence area; if they are close $\to$ small areas; if there are none $\to$ infinite area (but I've seen some polylogs so I doubt it will be the case :-(... This should help you estimate whether or not computing the series seems worth it.
• if it looks decent, compute the series, verify that if does converge where you want to; and then implement the method.

Hope this was useful.