8. Series solutions

Every technique so far assumed the equation had a nice enough shape: separable, linear, exact, or constant-coefficient. Many important equations have none of these shapes — including some of the most physically significant ones — and for those, Module 8's Taylor series machinery becomes a genuine solving technique in its own right, not just an approximation tool.

The method

Assume the solution can be written as a power series (§8.0) centered at x=0:

y=\sum_{n=0}^\infty a_nx^n

Differentiate term by term (§8.4 justified this is valid strictly inside the radius of convergence) to get y' and, if needed, y'' as series too, then substitute every series into the differential equation. Collecting like powers of x on both sides forces the coefficients to satisfy a recurrence relation — a formula expressing a_{n+1} (or some later term) in terms of earlier ones. Starting from a_0 (and a_1, for a second-order equation), the recurrence generates every coefficient in sequence, building up the solution's Taylor series one term at a time — exactly the reverse of §8.1's process, where the coefficients came from known derivatives of an already-known function.

Doing it in Python

Solving y'=xy by series, and confirming the recurrence reproduces the known exact solution from §13.1 (y=Ce^{x^2/2}, found there by separation):

import sympy as sp

x = sp.Symbol('x')
N = 8   # how many terms to generate

# recurrence: (n+1)a_{n+1} = a_{n-1}  =>  a_{n+1} = a_{n-1}/(n+1)
a = [sp.Integer(0)] * (N + 1)
a[0] = 1   # this plays the role of the arbitrary constant C
a[1] = 0   # forced by matching the x^0 coefficient
for n in range(1, N):
    a[n + 1] = sp.Rational(a[n - 1], n + 1)

print(f"coefficients a_0 through a_{N}: {a}")

series_solution = sum(a[n] * x**n for n in range(N + 1))
print(f"series solution: {series_solution}")

exact_series = sp.series(sp.exp(x**2 / 2), x, 0, N + 1).removeO()
print(f"Taylor series of e^(x^2/2): {exact_series}")
print(f"match: {sp.expand(series_solution - exact_series) == 0}")

Deriving the recurrence itself — substituting the series into y'=xy and matching coefficients, done symbolically to confirm the by-hand derivation:

import sympy as sp

x, n = sp.symbols('x n')
a = sp.IndexedBase('a')

# y = sum a_n x^n  =>  y' = sum (n+1) a_{n+1} x^n,  xy = sum a_{n-1} x^n (shifted)
# matching coefficients of x^n for n >= 1:  (n+1)*a[n+1] = a[n-1]
recurrence = sp.Eq((n + 1) * a[n + 1], a[n - 1])
print(f"recurrence relation: {recurrence}")
print("solving for a[n+1]:", sp.solve(recurrence, a[n+1])[0])

Worked example

Find the first several terms of the series solution to y'=xy, and confirm it matches the known exact solution.

Write y=\sum_{n=0}^\infty a_nx^n, so y'=\sum_{n=0}^\infty(n+1)a_{n+1}x^n (re-indexed so the power of x matches directly). The right side:

xy=x\sum_{n=0}^\infty a_nx^n=\sum_{n=0}^\infty a_nx^{n+1}=\sum_{n=1}^\infty a_{n-1}x^n

(re-indexed by shifting n\to n-1, so the sum starts at n=1 since the lowest power on this side is x^1).

Matching the x^0 coefficient (present only on the left side, since the right side starts at x^1): 1\cdot a_1=0\Rightarrow a_1=0.

Matching the x^n coefficient for n\ge1:

(n+1)a_{n+1}=a_{n-1}\ \Longrightarrow\ a_{n+1}=\frac{a_{n-1}}{n+1}

Since a_1=0, every odd-indexed coefficient is forced to 0 (a_3=\frac{a_1}3=0, a_5=\frac{a_3}5=0, and so on). Starting from a_0 (the free constant), the even-indexed coefficients build up:

a_2=\frac{a_0}2,\qquad a_4=\frac{a_2}4=\frac{a_0}8,\qquad a_6=\frac{a_4}6=\frac{a_0}{48}

y=a_0\left(1+\frac{x^2}2+\frac{x^4}8+\frac{x^6}{48}+\cdots\right)

\boxed{y=a_0\sum_{k=0}^\infty\frac{x^{2k}}{2^kk!}=a_0e^{x^2/2}}

Sanity check. This is exactly §13.1's separation-of-variables answer, y=Ce^{x^2/2} (with a_0 playing the role of C) — a genuinely independent method arriving at the identical family of solutions, confirming both techniques are correct. The pattern a_{2k}=\frac{a_0}{2^kk!} is recognizable as exactly the Maclaurin coefficients of e^{x^2/2} from §8.1 (substitute u=x^2/2 into e^u's series, \sum\frac{u^k}{k!}) — the series method rediscovered a known Taylor series purely from the recurrence, without ever being told what function it was building. ✓

Why this method matters beyond cross-checking

y'=xy was chosen because it's independently solvable, to verify the series method works. Its real value is for equations with no elementary solution at all. The classic example is Airy's equation, y''-xy=0 — despite its simple appearance, no combination of polynomials, exponentials, trig functions, and logarithms solves it. The identical series-substitution method still applies without modification, producing a valid recurrence and a genuine power series solution (the Airy function) — one of the standard special functions in physics (appearing in optics and quantum mechanics), defined by exactly this series, because no simpler closed form exists.

Your turn

1. Verify by direct substitution that a_{n+1}=\dfrac{a_n}{n+1} is the correct recurrence for y'=y (§13.0's own opening example), and confirm it reproduces y=a_0e^x.

2. Using the recurrence from problem 1, find a_4 in terms of a_0.

3. True or false: series solutions can only ever reproduce a function that already had an elementary closed form.

Solutions

1. y=\sum a_nx^n\Rightarrow y'=\sum(n+1)a_{n+1}x^n. Setting y'=y and matching the x^n coefficient: (n+1)a_{n+1}=a_n\Rightarrow a_{n+1}=\dfrac{a_n}{n+1}. Starting from a_0: a_1=a_0, a_2=\dfrac{a_1}2=\dfrac{a_0}2, a_3=\dfrac{a_2}3=\dfrac{a_0}6, and in general a_n=\dfrac{a_0}{n!} — exactly e^x's own Maclaurin coefficients from §8.1, so y=a_0\sum\dfrac{x^n}{n!}=a_0e^x, matching §13.0's stated solution y=Ce^x exactly.

2. a_4=\dfrac{a_0}{4!}=\boxed{\dfrac{a_0}{24}}.

3. False. This is precisely the point of the Airy's-equation remark: series solutions routinely produce genuinely new functions with no elementary closed form at all, and the resulting power series (subject to §8.0's radius-of-convergence analysis) is a completely legitimate way to define and work with such a function — computing values, derivatives, and integrals from the series directly, exactly as §8.3 used the e^x and \arctan x series to compute e and \pi despite never writing those constants down any other way.

Check yourself in code

Using the recurrence a_{n+1}=\dfrac{a_{n-1}}{n+1} (with a_0=1, a_1=0) for y'=xy, generate coefficients a_0 through a_8 and confirm the resulting series matches the Taylor series of e^{x^2/2}.

Print exactly this:

coefficients: [1, 0, 1/2, 0, 1/8, 0, 1/48, 0, 1/384]
match: True

import sympy as sp

x = sp.Symbol('x')
N = 8

a = [sp.Integer(0)] * (N + 1)
a[0] = 1
a[1] = 0
for n in range(1, N):
    a[n + 1] = sp.Rational(a[n - 1], n + 1)
print("coefficients: ...")

series_solution = sum(a[n] * x**n for n in range(N + 1))
exact_series = sp.series(sp.exp(x**2 / 2), x, 0, N + 1).removeO()
print("match: ...")

import sympy as sp

x = sp.Symbol('x')
N = 8

a = [sp.Integer(0)] * (N + 1)
a[0] = 1
a[1] = 0
for n in range(1, N):
    a[n + 1] = sp.Rational(a[n - 1], n + 1)
print(f"coefficients: {a}")

series_solution = sum(a[n] * x**n for n in range(N + 1))
exact_series = sp.series(sp.exp(x**2 / 2), x, 0, N + 1).removeO()
print(f"match: {sp.expand(series_solution - exact_series) == 0}")

A series solution assumes y=\sum a_nx^n, substitutes into the differential equation, and matches coefficients to build a recurrence relation that generates every a_n in sequence — verified on y'=xy against §13.1's exact answer, and applicable without modification to equations like Airy's y''-xy=0 that have no elementary solution at all, where the resulting series becomes the definition of a genuinely new special function.

Next: differential equations involving more than one unknown function at once — systems of equations, and the phase-plane picture that visualizes their solutions without solving them explicitly.