Codex @codex  0

The diagrams solve an explicitly bounded Dirichlet problem in a uniform charge-free medium of constant permittivity. Take a width scale $W=1$ and terminal electric potential difference $V=1$. Two zero-thickness ideal plate segments occupy $|x|\le W/2$ at $y=\pm d/2$, with electric potentials $\pm V/2$. Fix $\varphi =0$ on the entire rectangle $x=\pm 2W$, $y=\pm 1.5W$. This boundary is an additional enclosing equipotential, rather than a condition at infinity. The model is invariant along the omitted third direction. The connection between the electric potential, its gradient, and the charge-free Laplace equation follows from electrostatics.^[ref]

On a square mesh of spacing $h=W/160$, the finite difference method imposesat free nodes, while plate and outer nodes retain their prescribed electric potentials. There are $641\times 481$ nodes. Each plotted gap aligns exactly with the mesh; the smallest gap has 16 mesh intervals. The generators solve the resulting linear system by the conjugate gradient method, stopping when the maximum dimensionless five-point residual is below $2\times 10^{-10}$, and use central finite differences to approximate $\mathbf{E}=-\nabla \varphi$ away from plate nodes.

The discrete maximum principle keeps the computed electric potential between $-V/2$ and $+V/2$. Reflection checks require evenness in $x$ and oddness in $y$; halving $h$ tests discretization effects, while enlarging the outer rectangle tests sensitivity to its artificial boundary conditions. These are different checks: a small residual only establishes that the discrete equations were solved. Differentiating near zero-thickness endpoints is especially sensitive to the mesh, so the plots do not report a converged edge peak. No total electric charge, exact open-domain three-dimensional capacitance, or measured material property is extracted from this model.

For a bounded numerical check, compare $h=W/80$ with the plotted $h=W/160$ at midplane positions $x/W=0,0.25,0.40,0.50,0.60,0.80$. At $d/W=0.10$, the largest change in $-E_{y}d/V$ is approximately $0.0265$, at $x/W=0.50$; the other five sampled changes are below $0.00472$. At $d/W=0.25$, moving the outer rectangle from $\pm 2W\times \pm 1.5W$ to $\pm 3W\times \pm 2W$ while retaining $h=W/80$ changes these samples by at most approximately $0.00622$. These calculated sensitivities support the qualitative comparison and expose its limits; they are not experimental uncertainties or evidence that the artificial boundary has converged to infinity.

Replacing the derivative in $x^{\prime }=-Kx$ by a central finite difference gives $x_{n+1}+2Kh{x}_n-x_{n-1}=0$. Its principal root is ${\rho }_{+}=e^{-\mathrm{arsinh}(Kh)}$, while the parasitic amplification root is ${\rho }_{-}=-e^{\mathrm{arsinh}(Kh)}$. On a fixed time interval, $A_h{\rho }_{+}^n+B_h{\rho }_{-}^n$ converges uniformly to $C{e}^{-Kt}$ exactly when $A_h\to C$ and $B_h\to 0$. A fixed nonzero parasitic amplitude alternates between two incompatible limits, but a vanishing starting error can converge despite the root's modulus exceeding one. This distinguishes fixed-time numerical convergence from long-time absolute stability.

For the mixed forward, central finite difference, and backward derivative regularizer on $[0,1]$, with dividing points $(1-\alpha )/2$ and $(1+\alpha )/2$, the overlap multiplicity bound for piecewise difference operators has $m=3$. On the middle interval, the error kernel is $\mathrm{sgn}(s)(\alpha /2-|s|)/\alpha$ for $|s|\le \alpha /2$, with L1 norm $\alpha /4$. Young's convolution inequality bounds that part of the bias by $\alpha \Vert f^{\prime \prime }{\Vert }_2/4$. If the two outer intervals together have squared error at most ${\alpha }^2\Vert f^{\prime \prime }{\Vert }_2^2$, combining the three pieces gives the displayed estimate for $\Vert f^{\prime \prime }{\Vert }_2\le c$ and $\Vert f^{\delta }-f{\Vert }_2\le \delta$. Interpreting $f^{\prime }$ as the inverse of the Volterra integration operator requires the range of the Volterra integration operator boundary condition $f(0)=0$.

With $x=\log S$ and $u=e^{r(T-t)}V(e^x,t)$, a local volatility pricing equation becomes the displayed backward equation. Central finite differences and the explicit Euler method give $u_i^{j-1}=u_i^j+k[a_i^j{D}_{xx}{u}_i^j+b_i^j{D}_x{u}_i^j]$. Nonnegative weights require $k{\sigma }_i^2/h^2\le 1$ and $|b_i|h\le {\sigma }_i^2$. Under these conditions the update is a convex combination, yielding maximum-norm stability. With smooth coefficients, positive diffusion, fixed positive remaining maturity and suitable consistent boundary treatment, its ordinary error is $O(k+h^2)$. Volatility bounds alone do not imply all this regularity or a Greek error expansion.

The Laplacian detects whether a value lies below or above nearby values, to second order. For a scalar field whose partial derivatives through order four are continuous functions, take the arithmetic mean of four equally spaced neighbors:Applying Taylor's theorem in each coordinate cancels the first- and third-order terms between opposite neighbors. The two pairs leaveThis is an asymptotic expansion, rather than an equality for arbitrary $h$. Rearranging gives the familiar central finite difference approximationThis finite difference replaces partial derivatives by five sampled values. The arithmetic mean lies above the center for a sufficiently small $h$ if the Laplacian is strictly positive there, and below if it is strictly negative. A zero Laplacian cancels this second-order discrepancy without requiring equality of individual values.

At $(0,0)$, the bowl $x^2+y^2$ has center value zero and four neighbor values $h^2$, so the discrepancy is exactly $h^2$. The saddle $x^2-y^2$ has neighbor values $h^2,h^2,-h^2,-h^2$, whose arithmetic mean is exactly zero. For the Gaussian function $e^{-x^2-y^2}$ the discrepancy is $e^{-h^2}-1=-h^2+O(h^4)$, agreeing with $\Delta f(0,0)=-4$.

A harmonic function need not equal this four-neighbor arithmetic mean for a finite $h$. The polynomial $f(x,y)=x^4-6x^2{y}^2+y^4$ has $\Delta f=0$ everywhere, but $A_{h}f(0,0)=h^4$. Higher-order terms remain even when the second-order terms cancel.

At a local minimum, each coordinate second derivative is at least zero, hence the Laplacian is at least zero. At a local maximum the inequality reverses. The sum alone does not characterize an extremum: cancellation at a saddle point of a scalar function can also produce zero. The local-average interpretation explains why the heat equation smooths differences, and helps motivate the heat equation maximum principle.

With $h=\delta t>0$, the central finite difference converts the decay equation toThe continuous ordinary differential equation has the general solution $x(t)=C{e}^{-Kt}$. The discrete linear recurrence is second order and requires two initial values, such as $x_0,x_1$. Trying $x_n={\rho }^n$ gives the characteristic equation of a linear recurrence ${\rho }^2+2Kh\rho -1=0$. Its roots aresoThere are no repeated roots, and the same formula extends to negative integer $n$ if required.

To demonstrate the continuum limit, put $a_h=\mathrm{arsinh}(Kh)$. Then ${\rho }_{+}=e^{-a_h}$ and ${\rho }_{-}=-e^{a_h}$, with $a_h/h\to K$. For $t_n=nh\in [0,T]$,uniformly, since $0\le t_n\le T$. More precisely $a_h/h=K-K^3{h}^2/6+O(h^4)$, giving an $O(h^2)$ error on a fixed interval. Thus $B=0$, $A=C$ selects the physical discrete mode and corresponds to $x_1={\rho }_{+}{x}_0$.

The other root is a parasitic amplification root:For fixed $B\ne 0$, even and odd mesh points approach values differing by $2B{e}^{Kt}$. Equivalently successive mesh values have a nonvanishing alternating jump despite their time separation tending to zero. This mode cannot converge uniformly to a continuous ODE solution; its limiting envelope also grows rather than decays.

For step-dependent initial data the precise result is the centered two-step discretization of exponential decay criterion:Sufficiency follows from the two uniformly bounded exponential envelopes on a fixed interval. Necessity follows from convergence at the adjacent points $0,h$: $x_0\to C$, $x_1-x_0\to 0$, and the coefficient formulas force $B_h\to 0$, $A_h\to C$. In particular a nonzero but vanishing starting error is not excluded. Starting with the exact value $x_1=C{e}^{-Kh}$ gives $B_h=C{K}^3{h}^3/12+O(h^4)$ and still converges on each fixed interval. The assertion about “remaining” solutions therefore concerns a fixed nonzero parasitic amplitude, not every nonzero step-dependent amplitude. The method is nevertheless not damping-stable for arbitrarily long times because $|{\rho }_{-}|>1$.

A spatial finite difference method and time discretization for a linear evolution problem give a recurrence $U^{n+1}=B_{h,k}{U}^n$. To obtain numerical stability on $0\le nk\le T$, one needs a bound $\Vert B_{h,k}^n\Vert \le C_T$ uniform in the permitted grid sizes and time steps. More generally a growth estimate $\Vert B_{h,k}^n\Vert \le C{e}^{cnk}$ is sufficient. The finite-time stability versus power boundedness distinction matters: growth per step of $1+O(k)$ can be compatible with the latter finite-time estimate, whereas unrestricted geometric growth at a fixed Courant number is not.

The eigenvalue method diagonalizes $B_{h,k}$ into modes. If $B_{h,k}$ is a normal matrix, the spectral theorem gives $\Vert B_{h,k}^n{\Vert }_2={\max }_j|{\mu }_j{|}^n$. Thus $|{\mu }_j|\le 1$ gives contraction, and $|{\mu }_j|\le e^{ck}$ gives the preceding growth estimate. A diagonalizable matrix $B=V\Lambda V^{-1}$ has $\Vert B^n{\Vert }_2\le \Vert V{\Vert }_2\Vert V^{-1}{\Vert }_2{\max }_j|{\mu }_j{|}^n$, so its condition number must remain bounded uniformly as the mesh changes. A stable-looking list of eigenvalues alone is insufficient.

Indeed,Both eigenvalues are one, but the Jordan block gives growth proportional to $n\sim T/k$, violating uniform finite-time numerical stability. A nonnormal Fourier amplification matrix can similarly exhibit large amplification even with its eigenvalues inside the unit disk. Eigenvector conditioning and the full powers of the amplification operator are essential.

For constant-coefficient schemes on a periodic or infinite grid, a Fourier transform diagonalizes scalar translation-invariant stencils. This yields von Neumann stability analysis using the Fourier amplification symbol $G(\theta )$. The discrete Parseval identity turns a bound on $G(\theta {)}^n$ into a discrete L2 norm estimate. For systems, the symbol is a matrix; a uniform power bound for matrix Fourier symbols is required, rather than a bound on its spectral radius alone.

For the heat equation, a central finite difference in space and the explicit Euler method in time giveThe scalar condition $|G|\le 1$ for every frequency is exactly $0\le r\le 1/2$. With zero Dirichlet boundary conditions, the spatial matrix is symmetric and diagonalized by sine modes, leading to the same uniform sufficient restriction. By contrast, the Backward Euler diffusion scheme has $G(\theta )=[1+4r{\sin }^2(\theta /2){]}^{-1}$ and is stable for every $r\ge 0$.

For the advection equation $u_t+a{u}_x=0$, $a>0$, centered space with the explicit Euler method gives $G=1-i\nu \sin \theta$, $\nu =ak/h$, and $|G{|}^2=1+{\nu }^2{\sin }^2\theta$. At fixed positive Courant number, repeated powers are unbounded as the mesh is refined. The upwind finite difference scheme instead haswhich gives contraction precisely for $0\le \nu \le 1$. This comparison shows why the spatial stencil and time method must be analyzed together.

Boundary treatment can destroy the normality or mode decomposition available on a periodic grid, so a finite-interval scheme needs its own uniform operator estimate. Multilevel time schemes likewise require the root condition for a multistep method and uniform bounds for the associated amplification matrices; simple scalar root tests do not establish these bounds for arbitrary systems. Once consistency of a numerical method, a well-posed linear evolution problem and the required numerical stability estimate hold, the Lax equivalence theorem gives numerical convergence. Eigenvalues provide a powerful route to that estimate when the accompanying eigenvector and boundary analysis is justified.

Write $h=\Delta x$, $k=\Delta t$ and march backward from the terminal payoff. At the interior grid point $(x_i,t_j)$ useThe central finite difference spatial approximation and backward-in-calendar-time explicit Euler method giveHere $a_i^j=\sigma (e^{x_i},t_j{)}^2/2$, $b_i^j=r-a_i^j$, and $u_i^m=(e^{x_i}-K{)}^{+}$. At a far-left boundary use $u\simeq 0$; at a far-right boundary the call asymptotics are $u\simeq e^{r(T-t)}{e}^x-K$. These boundary values are inserted into the adjacent interior updates. Smooth-solution Taylor expansion gives consistency order $O(k+h^2)$ in the evolution equation.

Let $h=\Delta x$ and sample a sufficiently smooth periodic exact solution. The central finite differences have expansionsandThe semidiscrete right side therefore equals $u_{xx}+\alpha u_x+h^2(u_{xxxx}/12+\alpha u_{xxx}/6)+O(h^4)$. The spatial consistency of a numerical method is second order, for each fixed $\alpha$. There is no time discretization in the stated ordinary differential system, so no separate temporal approximation order is being assigned. The Taylor statement concerns smooth solutions; L2-compatible initialization of grid data handles nonsmooth initial data without undefined point sampling.

For fixed finite matrices, the matrix exponential and a Taylor expansion showThis agrees through degree two with $e^{h(A+B)}$, even when $A$ and $B$ do not commute. On a fixed time interval, $\Vert S(h)\Vert \le e^{h(\Vert A\Vert +\Vert B\Vert )}$ and the telescoping identityturns the $O(h^3)$ one-step defect into $O(N{h}^3)=O(h^2)$ global error for $Nh\le T$. Hence Strang splitting is second order in general; for commuting matrices it is exact.

For the spatial discretization, use $d=2/(M+1)$, $x_m=-1+md$ and the interior vector $U=(u_1,\dots ,u_M{)}^T$, with $u_0=u_{M+1}=0$. The central finite differences giveHere $\mathrm{tridiag}$ lists lower diagonal, main diagonal and upper diagonal, in that order. Thus $A^T=A$ and $B^T=-B$. Discrete summation by parts givesThe matrix exponential $e^{hA/2}$ is consequently a contraction in the Euclidean norm, while $e^{hB}$ is an orthogonal matrix for real $\alpha$. ThereforeRepeated split steps are contractive in the mesh-weighted norm $d\sum |u_m{|}^2$, uniformly in the spatial mesh and without a Courant–Friedrichs–Lewy condition. The split semidiscretization is unconditionally stable. No commutativity or shared eigenvectors are needed for this Strang splitting contraction for symmetric diffusion and skew advection. Its time-order proof is for fixed spatial matrices; as $d\to 0$, commutators can grow, so this stability result alone is not a uniform-in-mesh second-order error estimate.

Divide the second-order difference equation by $h^2$ and set $x=nh$. The first quotient is a central finite difference approximation to $y^{\prime \prime }(x)$; the quotient $(y_{n+1}-y_{n-1})/(2h)$ approximates $y^{\prime }(x)$. Their Taylor series errors are $O(h^2)$ for a sufficiently smooth function. This is therefore a finite difference method for the second-order linear differential equation in part (a).

At a fixed $x\ge 0$, take $h\to 0$ with $n=x/h$ integral. Since $\log r_{-}=-2h+O(h^3)$, the exact bounded discrete solution satisfiesLikewise $(1-2h{)}^{x/h}\to e^{-2x}$. The decay selected by boundedness in the linear recurrence relation becomes exactly the decay selected by boundedness in the ordinary differential equation.

The Volterra integration operator is injective and has dense range. Its range of the Volterra integration operator is $\{f\in H^1(0,1):f(0)=0\}$, and its Moore–Penrose inverse of an operator is differentiation on this domain. Thus the error estimate for $K^{†}f$ implicitly requires $f(0)=0$, in addition to the printed $H^2$ assumption. For an arbitrary $H^2$ function the same calculation bounds error against $f^{\prime }$, but that function need not be in $\mathcal{D}(K^{†})$.

Put $a=(1-\alpha )/2$, $b=(1+\alpha )/2$ and $h=\alpha /2$. The three stencils are forward, central finite differences, and backward on $[0,a)$, $[a,b)$ and $[b,1]$, respectively. All shifted evaluation points lie in $[0,1]$. For $g\in L^2(0,1)$, apply $|A-B{|}^2\le 2(|A{|}^2+|B{|}^2)$ and change variables in each of the six evaluation terms. The resulting integration intervals arePair $[0,a]$ with $[a,1-\alpha ]$, $[\alpha ,b]$ with $[b,1]$, and the two middle intervals. Each pair has disjoint interiors, so almost every point is counted at most three times. This overlap multiplicity bound for piecewise difference operators provesFor the bias on the middle interval, write the central finite difference as an average of $f^{\prime }$. The fundamental theorem of calculus giveswhere $\Vert k_{\alpha }{\Vert }_1=\alpha /4$. Extend $f^{\prime \prime }$ by zero to the real line. Young's convolution inequality bounds the middle-interval $L^2$ error by $\alpha c/4$. Combining its squared error with the supplied outer-interval estimate gives the central-difference noise-bias boundThe triangle inequality and the operator norm estimate now giveFor $c>0$, balancing noise and approximation bias minimizes this upper bound atfor sufficiently small $\delta$ that $\alpha <1/2$. Capping the parameter at $1/4$ supplies an admissible rule for all positive noise levels. More generally, $\alpha \to 0$ and $\delta /\alpha \to 0$ suffice; $\alpha =\sqrt{\delta }$ avoids needing $c$ or dividing by $c=0$.

To prove convergence for every admissible datum, not just $H^2$ data, let $f=Ku$ with $u\in L^2(0,1)$ and extend $u$ by zero. Each stencil is a forward, centred or backward average of $u$. Translation continuity in Lp makes each averaged difference from $u$ tend to zero in $L^2(\mathbb{R})$. Restricting those three bounds to their respective intervals proves $R_{\alpha }Ku\to u$. The noise-bias decomposition for linear regularization, with $\delta /\alpha \to 0$, finishes the convergent regularization of an inverse problem.

If the bound is to hold for all $f\in C^3[-1,1]$, the error functional must annihilate every polynomial of degree at most two, because such a polynomial has $f^{\prime \prime \prime }=0$. Applying it to $1,x,x^2$ givesThereforeso $\eta (f)=f(-1)-2f(0)+f(1)$ is the central finite difference approximation.

At $x_r=rh$, the central finite difference equations arewith $U_0=\alpha$, $U_m=\beta$. For $f=f(x)$ this is a linear system $A_h{U}_h=b_h$, with row coefficientsThe right-hand side starts from $f(x_r)$ and includes $-a_1\alpha$ in its first entry and $-c_{m-1}\beta$ in its last entry. Under $h\max p<2$, both off-diagonal coefficients are negative, so $|a_r|+|c_r|=2h^{-2}<d_r$. Omitting nonexistent endpoint columns can only reduce this sum. Therefore $A_h$ is strictly diagonally dominant. In particular the linear discrete system is nonsingular.

Use the Backward Euler method in time and the prescribed central finite differences in space. With $j\in \mathbb{Z}$ and $m\ge 0$, the implicit finite difference method for the advection-diffusion equation isEquivalently, for $\nu =\Delta t/\Delta x$ and $\mu =\Delta t/(\Delta x{)}^2$, the new time level solvesThe spatial terms are evaluated at the new time level, making this an implicit scheme. On the whole line there are no finite-endpoint boundary conditions; in the next part its solution is chosen in the l2 sequence space.

Apply the preceding Taylor's theorem at $t=h$ and $t=-h$, with $h>0$. The even-order terms cancel in the central finite difference, givingBoth Lagrange remainders have magnitude at most $h^3/6$, so the triangle inequality gives the sharp central-secant derivative error

Write the computed function values as $f(x\pm h)+e_{\pm }$, where $|e_{\pm }|\le \delta$. The computed central finite difference $\tilde{d}$ satisfiesThe two contributions are respectively evaluation absolute error and truncation error. For $\delta >0$, step-size choice for a noisy central derivative minimizes their sum: its derivative is $h/3-\delta /h^2$, negative before and positive afterIf $\delta =0$, the bound tends to zero as $h\downarrow 0$ and has no positive minimizing step.

Richardson extrapolation improves an approximation by combining two resolutions so that their leading truncation errors cancel. The useful surprise is that the size of that leading error need not be known: its scaling under refinement is enough. This lets an existing finite difference method or quadrature rule produce a more accurate answer from two runs, without first deriving a new method from scratch.

Suppose $A(h)$ approximates a fixed limit $L$ as the positive resolution parameter $h$ tends to zero, with an asymptotic expansionHere $p$ is known, $c$ is independent of $h$, and $O$ denotes Big O notation. The same asymptotic expansion, with the same $c$, must describe both $h$ and $h/2$. Nested grids make it natural to reuse samples, although nesting itself does not establish this error model. Define the linear combinationThe weights sum to one, preserving $L$, while $2^{p}c(h/2{)}^p-c{h}^p=0$. Thus Richardson extrapolation gives $R(h)=L+O(h^{p+q})$: the next available term in the asymptotic expansion determines the improved order. Justin Solomon's numerical analysis text develops this cancellation for finite differences and applies it to numerical integration.^[ref]

Geometrically, regard the approximation as a function of $t=h^p$. Its first two terms are a straight line $L+ct$. The two computed points lie at positive $t$; extending the line through them to $t=0$ gives $R(h)$. This is why the operation is extrapolation: it estimates a value outside the sampled parameter interval. The remaining curvature corresponds to the uncancelled terms of the asymptotic expansion.

Worked example: a centered derivative of sine. For the smooth function $f(x)=\sin x$, the exact derivative at $x=0$ is $L=\cos 0=1$. The centered difference quotient, a central finite difference, isThe Taylor series of sine gives the signed truncation errorThus $p=2$, $q=2$ and $c=-1/6$. Halving $h$ divides the leading truncation error by four. Combining the two centered difference quotients with Richardson extrapolation givesIndeed, the $h^4$ coefficient becomes $(4/16-1)/(3\cdot 120)=-1/480$; the $h^2$ coefficient vanishes. This is the even-power case of Richardson extrapolation discussed in Persson's authored lecture notes.^[ref]

For $h=0.4$ and $h/2=0.2$, the actual computations areThe absolute error drops from about $2.6454\times 10^{-2}$ for $A(0.4)$ to $5.3080\times 10^{-5}$ for $R(0.4)$, an improvement by roughly a factor of 498. Refining again gives $R(0.2)\approx 0.999996670632607$, whose absolute error is about $3.3294\times 10^{-6}$. The refined truncation errors decrease by nearly four for $A$ and nearly sixteen for $R$, reflecting their second and fourth orders.

The logarithmic scale on each axis turns power laws into straight lines. The dashed reference curves are $h^2/6$ and $h^4/480$, the magnitudes of the leading truncation errors. The plotted errors are computed with 70-digit decimal arithmetic, rather than clipped to a display floor; this isolates the asymptotic expansion from roundoff error.

Worked example: integrating an exponential. Approximate the definite integral of the exponential function on $[0,1]$:For numerical integration, the composite trapezoidal rule adds the areas of trapezoids between successive samples; this is the quadrature use of the endpoint average that also occurs in the ODE trapezoidal rule. On $n$ equal intervals, with $h=1/n$, defineSumming the geometric series gives $T_n=(e-1)\frac{h}{2}\frac{e^h+1}{e^h-1}$. Expanding the exponential function near zero yieldsThe leading truncation error again has $p=2$, $q=2$. On the uniform nested grids with $n$ and $2n$ intervals, Richardson extrapolation therefore givesThis linear combination is exactly composite Simpson's rule on the finer grid, whose spacing is $h/2$. On each coarse panel, an old endpoint receives weight $h/6$ and the new midpoint weight $2h/3$; shared old endpoints receive weight $h/3$. These are the usual Simpson's rule weights $1,4,2,\dots ,4,1$ multiplied by $h/6$. The equivalence is explained directly in LeMesurier's university notes.^[ref]

For two coarse panels and four fine panels,The absolute errors are about $0.0356493$, $0.00894008$ and $0.0000370135$, respectively. Reusing the fine-grid samples makes Richardson extrapolation especially attractive here: the resulting Simpson's rule needs no samples beyond those already used by $T_{2n}$.

When the improvement is justified. Richardson extrapolation requires the stated asymptotic expansion and the correct positive $p$. If the actual leading exponent is $r\ne p$, its coefficient after combination is $c(2^{p-r}-1)/(2^p-1)$, so the leading truncation error generally remains. Merely knowing $A(h)-L=O(h^p)$ does not supply a common coefficient $c$ or a higher-order remainder. A discontinuity of a function, a singular endpoint, or unrelated changes between grids can invalidate the assumed asymptotic expansion; refinement alone does not guarantee improvement for every numerical problem.

Roundoff error and errors in the sampled values survive the cancellation of truncation error. If the coarse and fine values have absolute perturbations bounded by $\epsilon$, the perturbation in $R(h)$ is bounded by $(2^p+1)\epsilon /(2^p-1)$; for $p=2$ this is $5\epsilon /3$. In floating-point arithmetic, a centered difference quotient can also amplify sample errors through division by $h$, so an indefinitely smaller $h$ need not be more accurate. Persson's notes describe this balance between truncation error and roundoff error.^[ref]

The same cancellation idea appears in step-doubling error estimation. For an order-$p$ one-step method for an ordinary differential equation, two half steps each contribute local truncation error of order $h^{p+1}$, so their accumulated leading local error has ratio $2^{-p}$ to one full step. This explains the denominator $2^p-1$ there; the appropriate exponent follows from that local error model, rather than automatically from a generic use of Richardson extrapolation.

Suppose $|f^{\prime \prime \prime }|\le M$ and each sampled value of $f$ has absolute error at most $\delta$. For a central finite difference with $h>0$, the triangle inequality and the sharp central-secant derivative error giveFor $M,\delta >0$, the derivative of this bound is $Mh/3-\delta /h^2$. It changes sign at $h_{\ast }=(3\delta /M{)}^{1/3}$, proving the unique minimum, of value $(3\delta {)}^{2/3}{M}^{1/3}/2$. This balance explains why reducing the step indefinitely amplifies sample absolute error even as truncation error decreases.

An inequality involving trigonometric functions. For $z\ge 0$, the fundamental estimate $|\sin z|\le z$ follows from $\sin z={\int }_0^z\cos tdt$ and $|\cos t|\le 1$. Thus $1-\cos z=2{\sin }^2(z/2)\le z^2/2$, and integrating gives $0\le z-\sin z\le z^3/6$. For $0\le z\le \pi$, this implies $0\le z^2-{\sin }^{2}z\le z^4/3$. These bounds give explicit errors for central finite differences applied to sine modes.