Unit D8: Series Solutions and Numerical MethodsD8 单元：级数解与数值方法

Suppose $y = \sum_{n=0}^{\infty} a_n x^n$. We wish to write $y'' $ as a series in $x^n$ so it can be added to $y$.

$$ y'' = \sum_{n=2}^{\infty} n(n-1)\, a_n x^{n-2}. $$

Let $k = n-2$, so $n = k+2$ and the lower limit $n=2$ becomes $k=0$:

$$ y'' = \sum_{k=0}^{\infty} (k+2)(k+1)\, a_{k+2}\, x^{k}. $$

Renaming the dummy index $k$ back to $n$, both $y''$ and $y$ are now series in $x^n$, so their sum has coefficient $(n+2)(n+1)a_{n+2} + a_n$ on $x^n$.

Suppose we must collect $y'' - 2x\,y'$ as a single series in $x^n$, where $y = \sum_{n=0}^{\infty} a_n x^n$. The two pieces start at different powers, so the bookkeeping matters.

From Example 1.1, $y'' = \sum_{n=0}^{\infty} (n+2)(n+1)a_{n+2}\, x^n$. For the second piece, $y' = \sum_{n=1}^{\infty} n\, a_n x^{n-1}$, so

$$ -2x\,y' = -2\sum_{n=1}^{\infty} n\, a_n x^{n} = \sum_{n=0}^{\infty} (-2n\, a_n)\, x^n, $$

where the lower limit may be written $n=0$ because the $n=0$ term contributes nothing. Adding,

$$ y'' - 2x\,y' = \sum_{n=0}^{\infty} \big[ (n+2)(n+1)a_{n+2} - 2n\, a_n \big] x^n. $$

The coefficient of $x^n$ is now in a single bracket, ready to be set equal to whatever appears on the right side of the equation. This is exactly the form needed for the Hermite-type equation $y'' - 2x y' + \lambda y = 0$.

Find the radius of convergence of $\sum_{n=0}^{\infty} \dfrac{x^n}{3^n (n+1)}$.

Ratio test. With $a_n = 1/[3^n(n+1)]$,

$$ \left| \frac{a_{n+1}}{a_n} \right| = \frac{3^n (n+1)}{3^{n+1}(n+2)} = \frac{1}{3}\cdot\frac{n+1}{n+2} \longrightarrow \frac{1}{3}, $$

so $1/R = 1/3$ and $R = 3$. The series converges for $|x| < 3$.

Sanity check by comparison. For $|x| < 3$ the factor $(x/3)^n$ decays geometrically and $1/(n+1) \le 1$, so the terms are dominated by a convergent geometric series; for $|x| > 3$ the geometric factor grows and the terms cannot tend to zero. Both routes agree on $R = 3$, which is the distance from the center $0$ to the singularity of the closed form $-\tfrac{3}{x}\ln(1 - x/3)$ at $x = 3$.

Take $x_0 = 0$, an ordinary point. Put $y = \sum_{n=0}^{\infty} a_n x^n$. Then $y'' = \sum_{n=0}^{\infty} (n+2)(n+1) a_{n+2} x^n$ after shifting. Substituting:

$$ \sum_{n=0}^{\infty} \big[ (n+2)(n+1) a_{n+2} + a_n \big] x^n = 0. $$

Each bracket must vanish, giving the recurrence

$$ a_{n+2} = -\frac{a_n}{(n+2)(n+1)}. $$

With $a_0, a_1$ free, the even coefficients reproduce $\cos x$ and the odd ones $\sin x$:

$$ y = a_0 \Big(1 - \tfrac{x^2}{2!} + \tfrac{x^4}{4!} - \cdots\Big) + a_1 \Big(x - \tfrac{x^3}{3!} + \cdots\Big) = a_0 \cos x + a_1 \sin x. $$

Airy's equation has no elementary solution, so the series is the answer rather than a detour. Take $x_0 = 0$, an ordinary point since the coefficients $p=0$, $q=-x$ are entire. Put $y = \sum_{n=0}^{\infty} a_n x^n$. Then

$$ y'' = \sum_{n=0}^{\infty} (n+2)(n+1)a_{n+2}\,x^n, \qquad x y = \sum_{n=1}^{\infty} a_{n-1} x^{n}. $$

Collecting, the constant term gives $2\cdot 1\cdot a_2 = 0$, so $a_2 = 0$. For $n \ge 1$,

$$ (n+2)(n+1)a_{n+2} - a_{n-1} = 0 \quad\Longrightarrow\quad a_{n+2} = \frac{a_{n-1}}{(n+2)(n+1)}. $$

This is a three-step recurrence, so the coefficients split into three chains by remainder of the index modulo $3$. Starting from $a_0$ (free) and $a_1$ (free), with $a_2 = 0$ forcing the entire $a_{2}, a_5, a_8, \dots$ chain to vanish:

$$ a_3 = \frac{a_0}{3\cdot 2} = \frac{a_0}{6}, \quad a_4 = \frac{a_1}{4\cdot 3} = \frac{a_1}{12}, \quad a_6 = \frac{a_3}{6\cdot 5} = \frac{a_0}{180}. $$

Hence the general solution is

$$ y = a_0\Big(1 + \tfrac{x^3}{6} + \tfrac{x^6}{180} + \cdots\Big) + a_1\Big(x + \tfrac{x^4}{12} + \cdots\Big), $$

the two independent series being the Airy functions. Because the coefficients are entire, the radius of convergence is infinite.

Solve $y'' - x y' - y = 0$ as a series about $x_0 = 1$, finding the recurrence. Substitute $t = x - 1$, so $x = t + 1$ and derivatives in $t$ equal derivatives in $x$. The equation becomes

$$ y'' - (t+1) y' - y = 0, $$

where primes are now $d/dt$. Put $y = \sum a_n t^n$. Then $y'' = \sum (n+2)(n+1)a_{n+2} t^n$, while

$$ -(t+1)y' = -\sum_{n} n\, a_n t^{n} - \sum_{n}(n+1)a_{n+1} t^{n}, \qquad -y = -\sum a_n t^n. $$

Collecting the coefficient of $t^n$ and setting it to zero,

$$ (n+2)(n+1)a_{n+2} - n\, a_n - (n+1)a_{n+1} - a_n = 0, $$ $$ a_{n+2} = \frac{(n+1)a_{n+1} + (n+1)a_n}{(n+2)(n+1)} = \frac{a_{n+1} + a_n}{n+2}. $$

With $a_0 = y(1)$ and $a_1 = y'(1)$ free, every later coefficient follows. The lesson is procedural: shifting the center to $x_0$ is identical to expanding the variable coefficients in powers of $t = x - x_0$ before collecting.

A two-term recurrence of the form $a_{n+2} = c_n\, a_n$ splits the coefficients into two independent chains: the even chain $a_0, a_2, a_4, \dots$ and the odd chain $a_1, a_3, a_5, \dots$. The even chain is fixed entirely once $a_0$ is chosen, and the odd chain once $a_1$ is chosen. Hence the solution space is spanned by the two series obtained from $(a_0,a_1)=(1,0)$ and $(0,1)$. These are linearly independent because at $x_0$ one has value $1$ and derivative $0$ while the other has value $0$ and derivative $1$, so their Wronskian at $x_0$ equals $1 \neq 0$. This matches the general theory: a second-order linear ODE has a two-dimensional solution space.

Solve $x^2 y'' - 3x y' + 4y = 0$ for $x>0$. Here $\alpha = -3$, $\beta = 4$, so the indicial equation is

$$ r(r-1) - 3r + 4 = r^2 - 4r + 4 = (r-2)^2 = 0, $$

a repeated root $r=2$. The general solution is therefore

$$ y = (c_1 + c_2 \ln x)\, x^2. $$

Solve $x^2 y'' + x y' + 4y = 0$ for $x > 0$ with $y(1) = 0$, $y'(1) = 2$.

Here $\alpha = 1$, $\beta = 4$, so the indicial equation is $r(r-1) + r + 4 = r^2 + 4 = 0$, giving $r = \pm 2i$, that is $\lambda = 0$, $\mu = 2$. The general solution is

$$ y = c_1 \cos(2\ln x) + c_2 \sin(2\ln x). $$

Apply the data at $x = 1$, where $\ln 1 = 0$: $y(1) = c_1 = 0$. Differentiating, $y' = -2c_1 \tfrac{\sin(2\ln x)}{x} + 2c_2 \tfrac{\cos(2\ln x)}{x}$, so $y'(1) = 2c_2 = 2$, giving $c_2 = 1$. Hence

$$ y = \sin(2\ln x). $$

The solution oscillates with respect to $\ln x$, so its zeros at $x = e^{k\pi/2}$ crowd together as $x \to 0^+$ and spread out as $x$ grows, the signature behavior of an Euler equation with complex exponents.

Solve $2x^2 y'' + 3x y' - y = 0$ for $x > 0$. Divide by $2$ to match the standard form $x^2 y'' + \tfrac{3}{2} x y' - \tfrac{1}{2} y = 0$, so $\alpha = \tfrac{3}{2}$, $\beta = -\tfrac{1}{2}$. The indicial equation is

$$ r(r-1) + \tfrac{3}{2} r - \tfrac{1}{2} = r^2 + \tfrac{1}{2} r - \tfrac{1}{2} = 0 \;\Longrightarrow\; 2r^2 + r - 1 = (2r - 1)(r + 1) = 0, $$

with roots $r_1 = \tfrac{1}{2}$ and $r_2 = -1$. The general solution is

$$ y = c_1 \sqrt{x} + c_2 x^{-1}. $$

The branch $x^{-1}$ blows up as $x \to 0^+$, so any solution bounded at the origin must have $c_2 = 0$. This is how Euler equations encode a boundary condition at a singular endpoint: only one of the two exponents is admissible there.

The substitution $t = \ln x$ (for $x>0$) turns an Euler equation into a constant-coefficient ODE. By the chain rule $\frac{dy}{dx} = \frac{1}{x}\frac{dy}{dt}$ and $\frac{d^2y}{dx^2} = \frac{1}{x^2}\big(\frac{d^2y}{dt^2} - \frac{dy}{dt}\big)$. Substituting into $x^2 y'' + \alpha x y' + \beta y = 0$ clears every power of $x$:

$$ \frac{d^2 y}{dt^2} + (\alpha - 1)\frac{dy}{dt} + \beta y = 0. $$

Its characteristic equation $r^2 + (\alpha-1)r + \beta = 0$ is identical to the indicial equation $r(r-1)+\alpha r + \beta = 0$. The three solution forms ($e^{rt}$, $t e^{rt}$, $e^{\lambda t}\cos\mu t$) become, on undoing $t=\ln x$, the three forms tabulated above. This explains why $\ln x$ appears in the repeated-root case: it is the image of $t$.

Approximate $y(0.2)$ for $y' = x + y$, $y(0)=1$, using $h=0.1$.

Step 1, from $(x_0,y_0)=(0,1)$: slope $f = 0 + 1 = 1$, so

$$ y_1 = 1 + 0.1\,(1) = 1.1, \qquad x_1 = 0.1. $$

Step 2, from $(0.1, 1.1)$: slope $f = 0.1 + 1.1 = 1.2$, so

$$ y_2 = 1.1 + 0.1\,(1.2) = 1.22, \qquad x_2 = 0.2. $$

So $y(0.2) \approx 1.22$. The exact solution is $y = 2e^x - x - 1$, giving $y(0.2) = 2e^{0.2} - 1.2 \approx 1.2428$, so the two-step estimate is low by about $0.023$.

Redo $y' = x + y$, $y(0) = 1$, this time reaching $x = 0.2$ in four steps of $h = 0.05$, and compare the error to the two-step run above.

So $y(0.2) \approx 1.23101$. The exact value is $1.24281$, an error of about $0.0118$. The two-step run ($h=0.1$) had error $0.0228$. Halving $h$ cut the error from $0.0228$ to $0.0118$, almost exactly the factor of two predicted for a first-order method, confirming the $O(h)$ scaling of Section 6.

The plain method uses only the slope at the left endpoint. Heun's method (the improved Euler, a second-order Runge-Kutta) averages the left slope with the slope at the Euler-predicted right endpoint:

$$ \tilde y_{n+1} = y_n + h f(x_n, y_n), \qquad y_{n+1} = y_n + \tfrac{h}{2}\big[f(x_n,y_n) + f(x_{n+1}, \tilde y_{n+1})\big]. $$

For $y' = x + y$, $y(0)=1$, $h = 0.2$: the left slope is $f(0,1) = 1$, the predictor is $\tilde y_1 = 1 + 0.2(1) = 1.2$, and the right slope is $f(0.2, 1.2) = 0.2 + 1.2 = 1.4$. Averaging,

$$ y_1 = 1 + \tfrac{0.2}{2}(1 + 1.4) = 1 + 0.1(2.4) = 1.24. $$

Against the exact $1.24281$ this single Heun step has error $0.0028$, an order of magnitude better than one plain Euler step of the same size. The extra slope evaluation buys an extra order of accuracy.

Take $y' = x + y$, $y(0)=1$, $h=0.2$, and compute $y_1 \approx y(0.2)$.

$$ k_1 = 0 + 1 = 1. $$ $$ k_2 = (0+0.1) + (1 + 0.1\cdot 1) = 0.1 + 1.1 = 1.2. $$ $$ k_3 = (0+0.1) + (1 + 0.1\cdot 1.2) = 0.1 + 1.12 = 1.22. $$ $$ k_4 = (0+0.2) + (1 + 0.2\cdot 1.22) = 0.2 + 1.244 = 1.444. $$ $$ y_1 = 1 + \tfrac{0.2}{6}\big(1 + 2(1.2) + 2(1.22) + 1.444\big) = 1 + \tfrac{0.2}{6}(7.284) = 1.2428. $$

This matches the exact $y(0.2)=1.24281\ldots$ to four decimals in a single step, far better than the two Euler steps of Example 4.1.

The second-order midpoint method is

$$ k_1 = f(x_n, y_n), \qquad k_2 = f\!\left(x_n + \tfrac{h}{2}, \, y_n + \tfrac{h}{2}k_1\right), \qquad y_{n+1} = y_n + h\, k_2. $$

Apply it to $y' = -2xy^2$, $y(0) = 1$, with $h = 0.2$, whose exact solution is $y = 1/(1+x^2)$. The left slope is $k_1 = -2(0)(1)^2 = 0$. The half-step predictor is $y_n + \tfrac{h}{2}k_1 = 1 + 0.1(0) = 1$ at $x = 0.1$, so

$$ k_2 = -2(0.1)(1)^2 = -0.2, \qquad y_1 = 1 + 0.2(-0.2) = 0.96. $$

The exact value is $1/(1 + 0.04) = 0.96154$, so the error is about $0.0015$. Plain Euler here would give $y_1 = 1 + 0.2(0) = 1$, missing the decrease entirely because the left slope is zero. The midpoint sample detects the curvature that Euler is blind to.

RK4 costs four slope evaluations per step. A fair comparison gives plain Euler the same budget: four Euler steps for one RK4 step. Take $y' = y$, $y(0) = 1$, target $y(0.4)$, exact value $e^{0.4} = 1.49182$.

One RK4 step, $h = 0.4$. With $f = y$: $k_1 = 1$, $k_2 = 1 + 0.2(1) = 1.2$, $k_3 = 1 + 0.2(1.2) = 1.24$, $k_4 = 1 + 0.4(1.24) = 1.496$. Then

$$ y_1 = 1 + \tfrac{0.4}{6}\big(1 + 2(1.2) + 2(1.24) + 1.496\big) = 1 + \tfrac{0.4}{6}(7.376) = 1.49173. $$

Four Euler steps, $h = 0.1$. Each step multiplies by $(1 + 0.1) = 1.1$, so $y_4 = 1.1^4 = 1.4641$.

RK4 errs by $0.00009$; Euler at equal cost errs by $0.0277$, roughly three hundred times worse. The weighted four-slope average is not just incrementally better, it is in a different accuracy class.

Suppose Euler's method with $h=0.1$ produces a global error of about $0.02$ at the endpoint. Estimate the error with $h=0.025$.

The step was reduced by a factor of $4$. For a first-order method the global error scales like $h^1$, so

$$ \text{new error} \approx 0.02 \times \frac{0.025}{0.1} = 0.02 \times \tfrac{1}{4} = 0.005. $$

For an RK4 run, the same factor-of-$4$ reduction would multiply the error by $(1/4)^4 = 1/256$, illustrating the payoff of higher order.

You run an unknown method twice at the same endpoint and record errors $E(h) = 4.0\times 10^{-3}$ at $h = 0.1$ and $E(h/2) = 5.0\times 10^{-4}$ at $h = 0.05$. What is the order $p$?

If $E(h) \approx K h^p$, then the ratio of errors over a halving is

$$ \frac{E(h)}{E(h/2)} \approx \frac{K h^p}{K (h/2)^p} = 2^p \quad\Longrightarrow\quad p = \log_2\!\frac{E(h)}{E(h/2)}. $$

Here the ratio is $4.0\times10^{-3} / 5.0\times10^{-4} = 8$, so $p = \log_2 8 = 3$. The method is third order. This log-ratio estimate of $p$ is the standard experimental check that a solver is implemented correctly: the measured order should match the theoretical one.

An RK4 integration over $[0,1]$ gives endpoint error $1.6\times 10^{-5}$ at $h = 0.1$. What step size $h^\star$ should reach a tolerance of $10^{-8}$?

For a fourth-order method $E \approx K h^4$, so $K \approx 1.6\times10^{-5} / (0.1)^4 = 0.16$. Solving $0.16 (h^\star)^4 = 10^{-8}$,

$$ (h^\star)^4 = \frac{10^{-8}}{0.16} = 6.25\times10^{-8}, \qquad h^\star = (6.25\times10^{-8})^{1/4} \approx 0.0158. $$

So roughly $h^\star \approx 0.016$, meaning about $1/0.016 \approx 63$ steps instead of $10$. Notice the strong leverage of high order: improving accuracy by a factor of $1600$ required shrinking $h$ by only about $6.3$, because $6.3^4 \approx 1600$.

Let $y(x)$ be the exact solution and expand it about $x_n$ by Taylor's theorem:

$$ y(x_n + h) = y(x_n) + h\, y'(x_n) + \frac{h^2}{2} y''(\xi), \quad \xi \in (x_n, x_n+h). $$

Since $y'(x_n) = f(x_n, y(x_n))$, the first two terms are exactly the Euler update $y_n + h f(x_n,y_n)$. The discrepancy in a single step is therefore

$$ y(x_n+h) - \big(y_n + h f(x_n,y_n)\big) = \frac{h^2}{2} y''(\xi) = O(h^2). $$

This is the local truncation error. Over an interval of fixed length there are about $1/h$ steps, so the errors accumulate to a global error of order $h^{2} \cdot h^{-1} = h^{1}$. That is precisely why Euler is called a first-order method, while its per-step error is second order.

Consider $(x^2 + 4) y'' + x y' + y = 0$ with a series sought about $x_0 = 0$. Writing it as $y'' + p y' + q y = 0$ gives $p = x/(x^2+4)$ and $q = 1/(x^2+4)$.

These are analytic except where $x^2 + 4 = 0$, that is at $x = \pm 2i$. The nearest such singularity to $x_0=0$ in the complex plane is a distance $2$ away. The theorem of Section 2 then guarantees the series solution converges at least for

$$ |x| < 2. $$

No computation of coefficients is needed for this guarantee; the radius is read directly from the location of the singular points.

Consider $2x y'' + y' + y = 0$, with a regular singular point at $x = 0$. Try the Frobenius series $y = \sum_{n=0}^{\infty} a_n x^{n+r}$, $a_0 \neq 0$. Then

$$ y' = \sum (n+r) a_n x^{n+r-1}, \qquad y'' = \sum (n+r)(n+r-1) a_n x^{n+r-2}. $$

Substituting, the lowest power $x^{r-1}$ comes only from $2xy''$ and $y'$:

$$ \big[2r(r-1) + r\big] a_0\, x^{r-1} = r(2r - 1) a_0\, x^{r-1}. $$

Since $a_0 \neq 0$, the indicial equation is $r(2r-1) = 0$, with roots $r = 0$ and $r = \tfrac{1}{2}$. The general coefficient relation, from the coefficient of $x^{n+r-1}$, is

$$ \big[2(n+r)(n+r-1) + (n+r)\big] a_n + a_{n-1} = 0 \;\Longrightarrow\; a_n = \frac{-a_{n-1}}{(n+r)(2n+2r-1)}. $$

The larger root $r = \tfrac{1}{2}$ always yields a genuine Frobenius solution; the smaller root $r = 0$ gives the second independent solution here because the roots do not differ by an integer. The Euler equation of Section 3 is the degenerate case where this recurrence forces $a_n = 0$ for $n \ge 1$, leaving the single term $x^r$.

Apply forward Euler to the test equation $y' = \lambda y$ with $\lambda < 0$, whose exact solution $e^{\lambda x}$ decays to zero. One step gives

$$ y_{n+1} = y_n + h \lambda y_n = (1 + h\lambda) y_n \;\Longrightarrow\; y_n = (1 + h\lambda)^n y_0. $$

The numerical solution decays only when the amplification factor satisfies $|1 + h\lambda| < 1$, that is

$$ -2 < h\lambda < 0 \;\Longrightarrow\; h < \frac{2}{|\lambda|}. $$

If $\lambda = -100$ (a fast-decaying component), forward Euler is forced to use $h < 0.02$ even though the true solution is nearly zero almost immediately. This is the essence of stiffness: stability, not accuracy, dictates a punishingly small step. Backward Euler gives $y_n = (1 - h\lambda)^{-n} y_0$, whose factor has magnitude below $1$ for every $h > 0$ when $\lambda < 0$, so it stays stable at any step size.

The order-drop seen for Euler is general, and the reason is the accumulation of local errors under a Lipschitz bound. Suppose a one-step method has local truncation error bounded by $\tau \le C h^{p+1}$, and that $f$ satisfies a Lipschitz condition $|f(x,u) - f(x,v)| \le L |u - v|$. Let $e_n = y(x_n) - y_n$ be the global error after $n$ steps. Subtracting the numerical update from the exact-solution update introduces one fresh local error and propagates the existing one through the method, giving the recursive bound

$$ |e_{n+1}| \le (1 + hL)\,|e_n| + C h^{p+1}. $$

Iterating this linear inequality from $e_0 = 0$ yields a geometric sum,

$$ |e_n| \le C h^{p+1} \sum_{j=0}^{n-1}(1 + hL)^{j} = C h^{p+1}\,\frac{(1+hL)^n - 1}{hL}. $$

Over a fixed interval of length $X = nh$ we have $(1 + hL)^n \le e^{nhL} = e^{LX}$, so

$$ |e_n| \le \frac{C\,h^{p}}{L}\big(e^{LX} - 1\big) = O(h^{p}). $$

One power of $h$ is lost to the factor $1/(hL)$ produced by summing about $1/h$ steps. This is why a method whose single-step error is $O(h^{p+1})$ converges globally at the rate $O(h^p)$, and why the constant grows with the interval length and the Lipschitz constant of $f$.