Practice: Numerical Solution of Differential Equations

(a) Four Euler steps. With $f(t,y)=y-t^2+1$, $h=0.5$, applying $y_{n+1}=y_n+hf(t_n,y_n)$: - $f(0,0.5)=0.5-0+1=1.5$, so $y_1=0.5+0.5(1.5)=1.25$
- $f(0.5,1.25)=1.25-0.25+1=2.0$, so $y_2=1.25+0.5(2.0)=2.25$
- $f(1,2.25)=2.25-1+1=2.25$, so $y_3=2.25+0.5(2.25)=3.375$
- $f(1.5,3.375)=3.375-2.25+1=2.125$, so $y_4=3.375+0.5(2.125)=4.4375$ | $n$ | $t_n$ | Euler $y_n$ | Exact | Error |
|---|---|---|---|---|
| 0 | $0$ | $0.500000$ | $0.500000$ | $0$ |
| 1 | $0.5$ | $1.250000$ | $1.425639$ | $0.175639$ |
| 2 | $1.0$ | $2.250000$ | $2.640859$ | $0.390859$ |
| 3 | $1.5$ | $3.375000$ | $4.009155$ | $0.634155$ |
| 4 | $2.0$ | $4.437500$ | $5.305472$ | $0.867972$ | The approximation to $y(2)$ is $4.4375$, low by $0.86797195$, about 16%. The error grows monotonically: each step adds a new error to those already carried, and since this solution curves upward, every tangent line falls below it. Check on the exact solution: $y(0)=1-0.5=0.5$, and $y^{\prime }=2(t+1)-\frac{1}{2}{e}^t$ while $y-t^2+1=(t+1{)}^2-\frac{1}{2}{e}^t-t^2+1=2t+2-\frac{1}{2}{e}^t$, identical. (b) One improved Euler step.

Predict the far endpoint with an Euler step: $y_0+h{k}_1=0.5+0.5(1.5)=1.25$. Evaluate the second slope there:

Average:

Comparison at $t=0.5$. Exact $1.425639$; Heun's error $0.050639$, Euler's $0.175639$, better by a factor of about $3.5$ for one extra evaluation of $f$. Where the slopes were taken. $k_1$ at the current point $(0,0.5)$; $k_2$ at $(0.5,1.25)$, the predicted endpoint. Evaluating $k_2$ at the old value instead would give $f(0.5,0.5)=1.25$ and the wrong answer $1.1875$, and would forfeit the second-order accuracy, since averaging the slope at the two ends is what cancels the leading error term.

(a) Order from the error ratio.

The ratio is near 2, so halving the step roughly halves the error: first order, $p=1$, as expected for Euler. Prediction at $n=16$. Dividing by about 2 gives roughly $0.076$ to $0.080$. The computed value is $0.08035333$, a ratio of $1.8978$. Why the ratios climb. Order is an asymptotic statement: the error behaves like $C{h}^p$ as $h\to 0$, with corrections of higher order in $h$ that are not negligible at coarse steps. Those corrections shrink as $h$ does, so the observed ratio approaches the limit from below: $1.53$, $1.69$, $1.82$, $1.90$, $1.95$, $1.97$, $1.99$ for $n$ up to 128. The same pattern appears for RK4, whose ratios run $10.63$, $13.01$, $14.42$, $15.19$ toward 16. A ratio short of the limit is evidence that $h$ is not yet small, not that the order is wrong. (b) Instability. Euler on $y^{\prime }=\lambda y$ gives $y_{n+1}=(1+h\lambda )y_n$, so after $n$ steps $y_n=(1+h\lambda {)}^n{y}_0$, a geometric sequence with ratio $1+h\lambda$. With $\lambda =-20$ and $h=0.2$ the ratio is $1-4=-3$, so the computed values are $1$, $-3$, $9$, $-27$, $81$, $-243$: growing because , and alternating in sign because the ratio is negative. The true solution $e^{-20t}$ decays monotonically and stays positive, reaching $2.061153622439\times {10}^{-9}$ at $t=1$. So this is not inaccuracy but qualitative failure: wrong sign, and wrong by eleven orders of magnitude. The arithmetic is correct; the recurrence has stopped tracking the equation. The condition. Boundedness requires

The largest usable step is anything strictly below $0.1$. The endpoint matters: at $h=0.1$ exactly the ratio is $-1$, giving $1,-1,1,-1,\dots$, which never decays, and the computed $y(1)$ is $+1.000000$. At $h=0.05$ the ratio is 0 and the method returns 0; at $h=0.025$, $9.09\times {10}^{-13}$; at $h=0.01$, $2.04\times {10}^{-10}$. (c) What caps the step for $y^{\prime }=-50y$. The stability limit, not the accuracy requirement. With $\lambda =-50$, boundedness needs . One decimal place is an undemanding accuracy target and would permit a much coarser step on a well-behaved equation, but at $h=0.2$ the ratio $1+h\lambda =-9$ makes the output grow and alternate in sign while the true solution decays to $e^{-50}\approx 2\times {10}^{-22}$. The bound depends on the decay rate alone. An equation with $\lambda =-1000$ forces however few digits are wanted. Such equations are called stiff, and they motivate implicit methods, which carry no such restriction. (d) Equal work, different methods. Comparing per step is misleading, because one RK4 step costs four evaluations of $f$ and one Euler step costs one. Evaluating $f$ is usually the expensive part, so the fair comparison fixes the number of evaluations, which is what the budget of 32 does. RK4 takes the coarsest step and still wins by a factor of about $8300$ over Euler and $340$ over improved Euler. The exponent beats the step size: error behaves like $C{h}^p$, so quadrupling $h$ costs a factor of $4^4=256$ for RK4 while the extra order buys far more. Doubling the budget widens the gap, since each doubling buys Euler a factor of 2, improved Euler about 3.9, and RK4 about 15.2. Circumstances that overturn the ranking. If $f$ is not smooth, the order theorem's hypotheses fail and RK4's four samples straddle the kink, collapsing its advantage toward first order. If the equation is stiff, the step is capped by stability rather than accuracy, the extra accuracy per step is wasted, and an implicit method beats all three. And for one or two digits by hand, Euler's simplicity can be worth more than RK4's accuracy. (e) Reading a computed answer. How many digits to trust. The coarsest two estimates differ by $1.9\times {10}^{-6}$, which measures the error of the coarser run rather than the finer. The evidence is the agreement across refinements: 100 steps gives $0.3502318445$, and 1000 and 10000 steps both give $0.3502318443$, unchanged in the digits shown. Agreement across three refinements supports ten decimal places: $y(1)=0.3502318443$. This is the standard technique when no exact solution exists, since there is nothing to compare against. What it provides. Values at chosen points, for one set of inputs, with an error that must be estimated rather than known. Here a ten-digit number for a quantity no formula in this corpus can express. What it does not provide. A function. Nothing here says how $y$ depends on $t$ in general, what happens for a different initial value, or how the answer shifts if a parameter changes; each requires another run. Nor does it reveal structure: a formula like $3e^{x^2}$ displays its growth rate and behaviour at infinity at a glance. A numerical method also never announces that the solution has ceased to exist, so for $y^{\prime }=y^2$ with $y(0)=1$, which blows up at $t=1$, the arithmetic returns numbers past that point just as readily, and they mean nothing.