# Is there something faster than AsymptoticDSolveValue for getting the terms in a power series solution?

Here's an example of a differential equation which Mathematica 13.1 just returns without solving

DSolve[{y''[x] + (5 x^3 + 3 x^2 + 5 x + 1) y'[x] + x y[x] == 0, y[0] == 1, y'[0] == 2}, y[x], x]


This specific differential equation isn't so important, it's just an example. If you want a power series solution around $$x = 0$$ you can do something like

AsymptoticDSolveValue[{y''[x] + (5 x^3 + 3 x^2 + 5 x + 1) y'[x] + x y[x] == 0, y[0] == 1, y'[0] == 2}, y[x], {x, 0, 1000}];


The problem is that I don't need 1,000 terms, I need more like 1,000,000 and that's going to take ages. I assume - perhaps wrongly - that if you painstakingly derived the recurrence relation for the coefficients in the power series solution, a command like RecurrenceTable[] would get the terms much faster, especially if we only need floats and not the exact values. But is there an automatic way to do this?

• I would go with finding the recurrence first. Dec 1, 2022 at 20:52

You can't obtain the recurrence formula with V 13.1

Looking at "Q&A with Calculus Developers: Live with the R&D team" Nov 30, 2022 at https://www.twitch.tv/videos/1666674721 for V 13.2, I did not see this feature being there (see around 7 minutes frame). Not even Maple can do this for series solutions. i.e. give explicit recurrence formula used for series solution.

I have a function myself that does this (but still in development) and it only handle linear 1st and second order odes'. For your ode, this is the recurrence formula it found

$$n=0$$ gives $$a_{2} = -\frac{a_{1}}{2}$$ $$n=1$$ gives $$a_{3} = -\frac{a_{0}}{6}-\frac{2 a_{1}}{3}$$ $$n=2$$ gives $$a_{4} = \frac{a_{0}}{24}+\frac{a_{1}}{4}$$ For $$3\le n$$, the recurrence equation is $$a_{n +2}= -\frac{5 n a_{n}}{\left(n +2\right) \left(n +1\right)}-\frac{\left(5 n -10\right) a_{n -2}}{\left(n +2\right) \left(n +1\right)}-\frac{\left(3 n -2\right) a_{n -1}}{\left(n +2\right) \left(n +1\right)}-\frac{a_{n +1}}{n +2}$$

Now you are able to generate very fast all terms you want. Here is the Mathematica code

a[2] = -a[1]/2
a[3] = -a[0]/6 - 2 a[1]/3
a[4] = a[0]/24 + a[1]/4
Table[a[n + 2] = -5*n*a[n]/((n + 2)*(n + 1)) - (5*n - 10)*
a[n - 2]/((n + 2)*(n + 1)) - (3*n - 2)*
a[n - 1]/((n + 2)*(n + 1)) - a[n + 1]/(n + 2), {n, 3, 10}];
sol = Sum[a[n]*x^n, {n, 0, 10}];
Collect[sol, {a[0], a[1]}]


Compare to

AsymptoticDSolveValue[{y''[x] + (5 x^3 + 3 x^2 + 5 x + 1) y'[x] +
x y[x] == 0}, y[x], {x, 0, 10}]


Note that I generates all the $$a_n$$ using normal Table command. You could try to use RecurrenceTable instead to see if performance is better.

You can now solve for $$a_0,a_1$$ easily given the initial conditions.

• Thanks, I was afraid of that but it is what it is. It's not TOO bad to substitute a[n] z^n into the differential equation and re-index appropriately, wish it was built-in though. Dec 1, 2022 at 23:43
• @Mr.G that will work for only ordinary points. For regular singular point you need to use Frobenius series, Now you can't just do what you said like this, as first you need to find the indicial equation. It gets complicated quickly. For irregular singular point, you have to use asymptotic methods. standard series methods no longer works. btw, the and re-index appropriately is not that easy. There are many edge cases to take care of and so on. Try to code it and you will see it is not trivial but still much easier for ordinary point that regular singular points. Your ode is ordinary point. Dec 2, 2022 at 0:07
• Oh sure, I completely believe the general case is a nightmare but luckily I don't have to deal with that. Dec 2, 2022 at 0:19
• The bottom of this section en.wikipedia.org/wiki/Symbolic_integration#Recent_advances seems to say that there is a Maple package algolib that can compute the recurrence formulas. That said it's just my vague understanding of the text and I could be wrong. Dec 17, 2022 at 8:24

You can obtain the coefficients using SeriesCoefficient on the DifferentialRoot of the differential equation

f = DifferentialRoot[
Function[{y,
x}, {y''[x] + (5 x^3 + 3 x^2 + 5 x + 1) y'[x] + x y[x] == 0,
y[0] == 1, y'[0] == 2}]]

rec = SeriesCoefficient[f[x], {x, 0, n}][[1, 1, 1, 0, 1]][y, n]


(* {5 n y[n] + (4 + 3 n) y[1 + n] + (10 + 5 n) y[2 + n] + (3 + n) y[ 3 + n] + (3 + n) (4 + n) y[4 + n] == 0, y[0] == 1, y[1] == 2, y[2] == -1, y[3] == -(3/2)} *)

In LaTeX :

$$\left\{5 n y(n)+(3 n+4) y(n+1)+(5 n+10) y(n+2)+(n+3) y(n+3)+(n+3) (n+4) y(n+4)=0,y(0)=1,y(1)=2,y(2)=-1,y(3)=-\frac{3}{2}\right\}$$

To compare with @Nasser 's recurrence we can use:

Solve[rec[[1]] /. n -> n - 2, y[2 + n]][[1, 1]] // Simplify


(* y[2 + n] -> -( 1/((1 + n) (2 + n)))(5 (-2 + n) y[-2 + n] + (-2 + 3 n) y[-1 + n] + 5 n y[n] + y[1 + n] + n y[1 + n]) *)

In LaTeX

$$y(n+2)\to -\frac{5 (n-2) y(n-2)+(3 n-2) y(n-1)+5 n y(n)+n y(n+1)+y(n+1)}{(n+1) (n+2)}$$