= Adaptive Euler step doubling for exponential decay
{c}
{title2=$y'=-y,\quad y(0)=1$}
Consider the nonstiff scalar <initial value problem> $y'=-y$, $y(0)=1$, on $0\leq t\leq8$. Its exact solution is the <exponential decay> $y(t)=e^{-t}$. Use <step-doubling error estimation> with the <explicit Euler method>: from the same numerical starting value $y_n$, form one coarse step $C_n$ and two half steps giving $F_n$:
$$
C_n=(1-h_n)y_n,\qquad F_n=(1-h_n/2)^2y_n,\qquad \widehat e_n=F_n-C_n=\frac{h_n^2y_n}{4}.
$$
Accept $y_{n+1}=F_n$ when the <absolute and relative error tolerances> test $E_n=|\widehat e_n|/s_n\leq1$ passes, with $s_n=10^{-5}+10^{-3}\max(|y_n|,|F_n|)$. There is no <Richardson extrapolation> in the accepted value. Both approximations have <order of a numerical method> one, while this <local error estimator> has step exponent $r=2$. Consequently the <adaptive step-size controller> uses a <square root>, not a cube root:
$$
q_n=\min\!\left(5,\max\!\left(0.2,0.9E_n^{-1/2}\right)\right),\qquad h_{\rm next}=q_nh_n.
$$
Start with $h_0=1$. A rejected trial retains both $t_n$ and $y_n$, caps $q_n$ at one, and retries. A zero estimate uses $q_n=5$. Before every trial, cap $h_n$ at $8-t_n$. Stop with a failure if <floating-point arithmetic> gives $t_n+h_n=t_n$, or after 100000 trials. These rules specify the example completely; they are not a prescription for every <ordinary differential equation> solver.
\Image[/adaptive-euler-decay-trajectory-and-steps.png]
{title=Adaptive two-half-step Euler decay: exact and numerical trajectories, accepted step sizes and rejected trials}
The first two trials, $h=1$ and $h=0.2$, have $E\approx247.525$ and $E\approx9.90099$ and are rejected. The next trial has $h\approx0.0572049$ and $E\approx0.81$ and is accepted. For positive decreasing numerical values, $s_n=10^{-5}+10^{-3}y_n$, so the <local error estimator> gives
$$
E_n=\frac{h_n^2}{4(10^{-3}+10^{-5}/y_n)}.
$$
When the relative term dominates, this suggests nearly constant accepted steps, approximately $2(0.9)\sqrt{10^{-3}}$. When $y_n$ falls below the crossover size $10^{-5}/10^{-3}=0.01$, the absolute term increasingly permits longer steps. The sharp final decrease in $h$ only lands on $t=8$; it does not indicate a sudden loss of smoothness.
The actual signed local defect from the same numerical input is
$$
\ell_n=y_ne^{-h_n}-F_n=y_n\left(\frac{h_n^2}{4}-\frac{h_n^3}{6}+O(h_n^4)\right).
$$
Here $O$ is <Big O notation>, obtained from <Taylor's theorem> for the <exponential function>. Thus $\widehat e_n$ agrees with the leading term as $h_n\to0$, but is not identical to $\ell_n$. At $h=1$, $y_n=1$, the estimate is $0.25$ while the actual local defect is $e^{-1}-1/4\approx0.117879$. Their closeness improves for the smaller accepted steps. Had we added $\widehat e_n$ using <Richardson extrapolation>, the accepted formula would instead be $y_n(1-h_n+h_n^2/2)$, of <order of a numerical method> two; the same raw difference would still scale as $h_n^2$ and would not estimate that improved value's leading $h_n^3$ defect.
\Image[/adaptive-euler-local-estimates-and-global-errors.png]
{title=Estimated versus actual local Euler defects, rejected trials, and actual accumulated global error scaled by the local tolerance}
The <global error> $g_n=y_n-e^{-t_n}$ obeys an exact <recurrence relation> involving the actual local defects:
$$
g_{n+1}=e^{-h_n}g_n-\ell_n,\qquad g_0=0.
$$
Each previous <global error> contracts, yet the local defects keep adding. In the plotted run, 101 steps are accepted and two trials are rejected. Every accepted $E_n$ is at most approximately $0.81$, while the largest $|g_{n+1}|/s_n$ is approximately $38.6927$. At $t=8$, the numerical value is approximately $0.00027184218$, compared with the exact $e^{-8}\approx0.00033546263$, giving <absolute error> approximately $6.36204\times10^{-5}$. The right-hand plot deliberately scales actual <global error> by the same $s_n$ used for local acceptance; it is not a new <global error> tolerance test.
Tightening both <absolute and relative error tolerances> by a factor ten gives 315 accepted steps and endpoint <absolute error> approximately $2.14922\times10^{-5}$; tightening both by a factor one hundred gives 991 accepted steps and endpoint <absolute error> approximately $6.94171\times10^{-6}$. This is a convergence comparison for this particular <initial value problem>. The endpoint error decreases by about a factor three per tenfold tightening, rather than by ten: a first-order <Euler method> accumulates defects across many steps, while its <local error estimator> varies quadratically with $h$. A <variable-step global ODE error bound> needs propagation information as well as accepted local tests.
The <absolute stability> boundary remains a separate constraint. On the <Dahlquist test equation> $y'=-\lambda y$ with $\lambda>0$, the accepted two-half-step formula has <stability function> $R(-\lambda h)=(1-\lambda h/2)^2$. It strictly damps perturbations only for $0<\lambda h<4$. Its <local error estimator> is $\lambda^2h^2y_n/4$, which can pass an absolute-tolerance test at tiny $y_n$ even outside this damping interval. In a <stiff two-mode linear system>, a fast mode can therefore constrain the step long after its visible amplitude has become small. An <adaptive step-size controller> does not turn an <explicit Euler method> into an <A-stable> method. Official solver documentation likewise distinguishes explicit <Runge Kutta methods> for nonstiff problems from implicit <Runge Kutta methods> and <backward differentiation formulas> for <stiff equations>.<SciPy 1.16.2 documentation/scipy.integrate.solve_ivp>[solve_ivp, method parameter: explicit versus implicit methods and stiffness]{ref}
Back to article page
