# Recursive computation of symbolic series

I don't have experience in Mathematica and I would like to use the ability it has to do symbolic computations to find the terms of a series of functions recursively. The terms I want to find are given by:

$$S_{n+1}(x,t)=i\int_0^{t}\left(\frac{\partial^2S_n(x,\tau)}{\partial{x^2}}-A_n(x,\tau)\right)d\tau$$ $$A_n(x,t)=\sum_{i=0}^nS(x,t)_{i}S(x,t)_{n-i}$$

I have written the following:

S[0, x_, t_] = Exp[-(x)^2]
A[0, x_, t_] = S[0, x, t]*S[0, x, t]

S[n_Integer, r_, t_] :=
S[n, r, t] =
Refine[Block[{e},
I*Integrate[
1/(2) D[S[n - 1, r, e], {r, 2}] - A[n - 1, r, e], {e, 0, t}]],
Assumptions ->
Element[r | t | \[Sigma] | l | \[Alpha] | \[Beta], Reals]];

A[n_Integer, y_, t_] :=
A[n, y, t] =
Refine[Sum[S[i, y, t]*S[n - i, y, t], {i, 0, n}],
Assumptions ->
Element[y | t | \[Sigma] | l | \[Alpha] | \[Beta], Reals]];

SN[n_Integer, r_, t_] := Sum[S[i, r, t], {i, 0, n}]

h = SN[5, r, t]


However, when I use this type of construction for more complicated equations or when I want to get higher-order terms, the running time starts to be of the order of hours. Is there a way to make this code more optimized?

• Seems you are using the same label "S" to define two functions (the first one and the one that depends on n_integer. Also, the first "S" function do not depend on the variable "t". Jul 2, 2021 at 14:16
• I would recommend to label different functions with different names, e.g. S and Sn, A and An, to avoid recurring errors/calls to the same function. Jul 2, 2021 at 14:18
• The first one is like "the initial condition" so that must be labeled with n=0 in order to sum all the terms at the end. I don't know if that is what you were referring to. Jul 2, 2021 at 14:51

The trick is to use partial memoization: only memoize the index $$n$$ but not the function parameters $$(x,t)$$,

Clear[S, A, SN];
S = Function[{x, t}, Exp[-x^2]];
A[n_Integer?NonNegative] := A[n] = Function[{x, t},
Evaluate@Sum[S[i][x,t] * S[n-i][x,t], {i, 0, n}]]
S[n_Integer?Positive] := S[n] = Function[{x, t},
Evaluate[Expand[I*Integrate[D[S[n-1][x,τ], {x,2}] - A[n-1][x,τ], {τ, 0, t}]]]]
SN[n_Integer?NonNegative] := SN[n] = Function[{x, t},
Evaluate@Sum[S[i][x,t], {i, 0, n}]]


Now evaluation is pretty quick:

SN[x, t]
(*    E^-x^2 - I E^(-2 x^2) t - 2 I E^-x^2 t - E^(-3 x^2) t^2 - 4 E^(-2 x^2) t^2 -
6 E^-x^2 t^2 + I E^(-4 x^2) t^3 + 6 I E^(-3 x^2) t^3 + 56/3 I E^(-2 x^2) t^3 +
20 I E^-x^2 t^3 + E^(-5 x^2) t^4 + 8 E^(-4 x^2) t^4 + 110/3 E^(-3 x^2) t^4 +
96 E^(-2 x^2) t^4 + 70 E^-x^2 t^4 - I E^(-6 x^2) t^5 - 10 I E^(-5 x^2) t^5 -
892/15 I E^(-4 x^2) t^5 - 3628/15 I E^(-3 x^2) t^5 - 2656/5 I E^(-2 x^2) t^5 -
252 I E^-x^2 t^5 + 4 I E^-x^2 t x^2 + 12 E^(-2 x^2) t^2 x^2 + 24 E^-x^2 t^2 x^2 -
68/3 I E^(-3 x^2) t^3 x^2 - 368/3 I E^(-2 x^2) t^3 x^2 - 120 I E^-x^2 t^3 x^2 -
106/3 E^(-4 x^2) t^4 x^2 - 952/3 E^(-3 x^2) t^4 x^2 - 1008 E^(-2 x^2) t^4 x^2 -
560 E^-x^2 t^4 x^2 + 248/5 I E^(-5 x^2) t^5 x^2 + 3088/5 I E^(-4 x^2) t^5 x^2 +
17304/5 I E^(-3 x^2) t^5 x^2 + 38848/5 I E^(-2 x^2) t^5 x^2 + 2520 I E^-x^2 t^5 x^2 -
8 E^-x^2 t^2 x^4 + 224/3 I E^(-2 x^2) t^3 x^4 + 80 I E^-x^2 t^3 x^4 +
808/3 E^(-3 x^2) t^4 x^4 + 3776/3 E^(-2 x^2) t^4 x^4 + 560 E^-x^2 t^4 x^4 -
9872/15 I E^(-4 x^2) t^5 x^4 - 91696/15 I E^(-3 x^2) t^5 x^4 - 74112/5 I E^(-2 x^2) t^5 x^4 -
3360 I E^-x^2 t^5 x^4 - 32/3 I E^-x^2 t^3 x^6 - 320 E^(-2 x^2) t^4 x^6 -
448/3 E^-x^2 t^4 x^6 + 6688/3 I E^(-3 x^2) t^5 x^6 + 38144/5 I E^(-2 x^2) t^5 x^6 +
1344 I E^-x^2 t^5 x^6 + 32/3 E^-x^2 t^4 x^8 - 15872/15 I E^(-2 x^2) t^5 x^8 -
192 I E^-x^2 t^5 x^8 + 128/15 I E^-x^2 t^5 x^10    *)


If you prefer an interface of the form S[n, x, t] etc., you can define these:

S[n_, x_, t_] := S[n][x, t]
A[n_, x_, t_] := A[n][x, t]
SN[n_, x_, t_] := SN[n][x, t]