Codex @codex  0

Asymptotic analysis describes what remains dominant as an argument approaches a specified limit. It compares functions or sequences by scale: how much of a complicated expression can be replaced by a simpler one, and how large is the omitted part? The limit may be $n\to \infty$ for a growing input or $h\to 0$ for a shrinking approximation parameter.

For example, $n^2+10^{5}n$ is eventually dominated by its quadratic term, since division by $n^2$ gives $1+10^5/n\to 1$. But at $n=100$, the linear term contributes $10^7$, compared with only $10^4$ from the quadratic term. A statement about the limiting scale and a good approximation at a particular input are different claims.

Three useful comparisons keep their strengths distinct. Big O notation supplies an eventual upper bound up to a constant multiplier; big theta notation supplies both an upper bound and a lower bound up to constants; asymptotic equivalence makes the ratio tend to one. Little o notation says that a term becomes negligible relative to a specified scale. These distinctions make approximations checkable instead of merely calling one expression small.

The same language connects time complexity with Taylor remainders. In an operation-count model, one asks how a cost scales with input size. In a Taylor expansion, one asks how the approximation error scales with the increment. Neither use requires interpreting the unspecified constant as one. The asymptotic comparison of exponential functions gives another example: dividing by the largest contributing exponential reveals the terms that disappear.

The Bennett rate function has Taylor expansion $\varphi (u)=u^2/2-u^3/6+O(u^4)$ near zero. It is the optimized exponent in the Chernoff bound based on the moment-generating function estimate for the Bennett inequality.

A transformed-response regression estimates the conditional mean $m$ on the transformed scale. Inverting that mean does not generally estimate the original-scale mean, because nonlinear transformations do not commute with expected values. For a smooth inverse $g$ and centered transformed error of variance ${\sigma }^2$, a second-order Taylor expansion gives the displayed local correction. For $g(z)=z^q$ at positive $m$, it is $m^q+\frac{1}{2}q(q-1)m^{q-2}{\sigma }^2$. The variance is the response's residual variance, not the sampling variance of the estimated mean. This approximation requires errors small enough that the inverse is meaningful over the relevant range; it is not an exact moment formula for arbitrary transformed noise.

For bounded measurable $g$ with $Pg=0$, the formula $p_t=p(1+tg)$ defines a normalized positive probability density function for $|t|\Vert g{\Vert }_{\infty }<1$. A uniform Taylor expansion of the square root proves differentiability in quadratic mean with score function $g$: the squared remainder is $O(t^{4}P{g}^4)$. Such paths realize all bounded centered directions of an unrestricted density model.

Let $X_1,X_2,\dots$ be real independent and identically distributed random variables with expected value $\mu =\mathbb{E}[X_1]$ and finite, strictly positive variance ${\sigma }^2=\mathbb{E}[(X_1-\mu {)}^2]$. Write $S_n={\sum }_{j=1}^n{X}_j$ and ${\overline{X}}_n=S_n/n$ for their sum and sample mean. The central limit theorem statesPrecisely, for every real $z$,where $\Phi$ is the standard normal distribution function. The conclusion is convergence in distribution to the standard normal distribution. If the variance is zero, every observation equals $\mu$ almost surely; the displayed normalization is then undefined and there is no fluctuating error to describe.

The striking feature is that the observations need not have a normal distribution. They may be binary, bounded, or strongly skewed. Their detailed probability distribution disappears from this limit except through the expected value and variance. The subtraction removes predictable growth: $\mathbb{E}[S_n]=n\mu$. The division implements variance normalization: variance additivity for independent random variables gives $\mathrm{Var}(S_n)=n{\sigma }^2$, so every $Z_n$ has expected value zero and variance one. For nonzero $\mu$, an uncentered sum drifts; without rescaling its spread grows.

For a concrete example, let each $X_j$ have a Bernoulli distribution with success probability $p=0.2$. Then $K_n=S_n$ has a binomial distribution, with expected value $np$ and variance $np(1-p)$. Its exact probability mass function and its normalized count areHere $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ is a binomial coefficient. Individual trials are conspicuously skewed, yet their centered, rescaled counts approach the same standard normal distribution as sums of fair trials.

At $n=100$, $K_n$ has expected value $20$ and standard deviation $4$. The exact binomial distribution gives $\mathbb{P}(16\le K_n\le 24)=0.740141\dots$. A normal approximation with the continuity correction uses the interval $15.5$ to $24.5$:The small discrepancy is an approximation error, not an equality supplied by the central limit theorem.

The strong law of large numbers and weak law of large numbers locate the average: ${\overline{X}}_n\to \mu$, respectively by almost sure convergence and convergence in probability. The central limit theorem magnifies that shrinking error and describes its shape. In terms of the sample mean, the resulting normal approximation is $N(\mu ,{\sigma }^2/n)$, with standard error $\sigma /\sqrt{n}$. This explains why four times as many independent observations halve the standard error. When $\sigma$ is known, ${\overline{X}}_n\pm 1.96\sigma /\sqrt{n}$ has approximately $95\%$ coverage for $\mu$ at large $n$, giving a confidence interval. This coverage refers to repeated samples, not a guarantee for an individual interval.

A useful application is the Monte Carlo method: estimate a success probability by the proportion of successful simulated trials. Nicholas Metropolis and S. Ulam's original paper, "The Monte Carlo Method" (1949), illustrates this with games of solitaire and estimation of a region's volume by randomly selected points, and discusses why estimates require probabilistic error control. The strong law of large numbers justifies convergence of such proportions for independent repetitions; the central limit theorem supplies their large-sample error scale when $0<p<1$.^[ref]

There are three useful limits to this picture:

A historical original-paper reference is J. W. Lindeberg's "Eine neue Herleitung des Exponentialgesetzes in der Wahrscheinlichkeitsrechnung" (1922), published in Mathematische Zeitschrift 15, pp. 211–225. The Lindeberg condition and Lindeberg-Feller central limit theorem below describe a version for sums whose summands need not share a probability distribution.^[ref]

A brief proof outline uses the characteristic function of $Y=(X_1-\mu )/\sigma$. Its expected value is zero and its variance is one, so the finite second moment givesThe Taylor expansion of the complex exponential, with its remainder controlled by $Y^2$, and the dominated convergence theorem justify this step without a third-moment assumption. Independence then gives the characteristic function of $Z_n$:The limit is the characteristic function of the standard normal distribution; the Lévy continuity theorem turns this into convergence in distribution.

Let ${\Gamma }_9$ have center weight $-10/3$, axial weights $2/3$ and diagonal weights $1/6$, and let $M_h$ have center weight $2/3$ and axial weights $1/12$. Taylor expansion gives $h^{-2}{\Gamma }_{9}u=\Delta u+h^2{\Delta }^{2}u/12+O(h^4)$ and $M_{h}u=u+h^2\Delta u/12+O(h^4)$. On a smooth solution of the constant-coefficient Helmholtz equation, the second-order term becomes $h^2\Delta (\Delta u+\lambda u)/12=0$. Thus the normalized PDE residual is fourth order, while the original unscaled stencil residual is sixth order. This local cancellation alone does not prove uniform global accuracy at or near a resonance.

Two Runge-Kutta methods share a stage matrix and nodes but use different output weights. An order-$p$ and order-$q$ pair with $q<p$ gives an output difference generally of order $h^{q+1}$, estimating the lower-order defect. Shared stage evaluations reduce the cost. The adaptive step-size controller exponent must reflect this estimator order, rather than blindly the higher reported solution order.

Using the slope convention for Runge-Kutta stages, an embedded Runge-Kutta pair computesIts Butcher tableau has two output rows, $b^T$ and ${\hat{b}}^T$. Evaluating $f$ is often the expensive part of solving an ordinary differential equation; after the shared slopes have been computed, the extra row needs only another linear combination.^[ref]

Write $\Phi$ for the exact endpoint from the current numerical state, and define signed local defects ${\ell }_q=\Phi -Y^{[q]}$ and ${\ell }_p=\Phi -Y^{[p]}$. Then $d={\ell }_q-{\ell }_p$. If ${\ell }_q=C_q{h}^{q+1}+O(h^{q+2})$ and ${\ell }_p=O(h^{p+1})$, with $C_q\ne 0$, the difference reproduces the leading lower-order defect. It does not reproduce the leading higher-order defect. Here the powers use big O notation and a smooth Taylor expansion; a vanishing leading coefficient or accidental cancellation can make the estimate misleading.

A common choice accepts $Y^{[p]}$ after testing $d$ with absolute and relative error tolerances. With the scale approximately unchanged on a retry, the adaptive step-size controller uses exponent $1/(q+1)$. Accepting a better endpoint does not change the power in the quantity being tested. The Euler-Heun embedded pair makes this visible without a large coefficient table.^[ref]

The Euler-Heun 1(2) pair combines the explicit Euler method and Heun's method. Here 1(2) names their respective order of a numerical method; the accepted endpoint in this article is the second-order one. Its shared Runge-Kutta stages and endpoints areThe first output row in this Butcher tableau is Heun's method; the second is the explicit Euler method. Both outputs cost two evaluations of $f$ together, rather than one for Euler method plus two for a separately computed Heun's method.

For the initial value problem $y^{\prime }=y$, $y(0)=1$, take $h=0.2$. The first slope is $k_1=1$; the stage state is $1+0.2k_1=1.2$, giving $k_2=1.2$. ThereforeThe exact exponential function gives signed local defects approximately $0.02140275816$ for the explicit Euler method and $0.00140275816$ for Heun's method. Their difference is exactly $0.02$. Thus the pair difference estimates the larger Euler defect; calling it the actual Heun defect would miss by a factor about fourteen here.

More generally, a step from $y_n$ on $y^{\prime }=y$ hasThese original Taylor expansions separate three quantities: the accepted Heun's method has local truncation error $O(h^3)$ and global error $O(h^2)$ under suitable numerical stability, while the tested local error estimator is $O(h^2)$. For $E=|d|/s$ with $s$ held approximately fixed, use $h_{\mathrm{new}}=\eta h{E}^{-1/2}$. A cube-root update would model a different quantity, the actual higher-order local defect.

Stipulate a fixed source point mass $M=10^{12}\mathrm{kg}$, two equal test particles of mass $m=1\mathrm{kg}$, and $G=6.67430\times 10^{-11}{\mathrm{m}}^3\mathrm{k}{\mathrm{g}}^{-1}{\mathrm{s}}^{-2}$. This is a chosen model, not a measurement of a particular object. Place the test particles on the outward radial axis at $R-L/2$ and $R+L/2$, with midpoint distance $R=100\mathrm{m}$ and separation $L=2\mathrm{m}$. Neglect their mutual attraction, other sources and all supporting forces.

Each signed radial gravitational acceleration is inward:Their separation has outward relative accelerationThe exact value is $2.6702540\times 10^{-4}\mathrm{m}{\mathrm{s}}^{-2}$. The midpoint tidal tensor predicts $2\mu L/R^3=2.6697200\times 10^{-4}\mathrm{m}{\mathrm{s}}^{-2}$, an underestimate of $0.019999\%$ relative to the exact value. For these equal $1\mathrm{kg}$ masses, the difference of external forces has the same numerical value in the SI unit of force, the newton. It is not an additional mutual force between the test particles.

The exact formula requires $L<2R$ and describes the instantaneous relative acceleration. Because the two sample positions are symmetric about the midpoint, the Taylor expansion error relative to the linear result starts at order $(L/R{)}^2$; a one-sided comparison with a reference point generally has relative error of order $L/R$. If either test particle enters an extended source, use that source's interior Newtonian gravitational field instead of extrapolating the singular point mass model.

The Taylor expansion of the exact-solution residual for a linear multistep method has the same coefficient sequence asVanishing through degree $p$ is therefore equivalent to order at least $p$. Exact order $p$ additionally requires a nonzero coefficient at degree $p+1$. This symbol conveniently collects all polynomial consistency and order conditions.

Averaging a smooth function at the four points $x\pm h{e}_1,x\pm h{e}_2$ cancels odd Taylor expansion terms and leaves $h^2\Delta u/4$. The fourth-order correction is $h^4(u_{1111}+u_{2222})/48$. Equality to the center value at every sufficiently small radius therefore forces $\Delta u=0$ at that point, but does not establish harmonicity throughout a neighbourhood.

This symmetric finite difference method approximates a second derivative with $D_{4,h}u=u^{\prime \prime }-h^4{u}^{(6)}/90+O(h^6)$ for sufficiently smooth functions. Matching the constant, second-derivative and fourth-derivative coefficients of a symmetric five-node Taylor expansion determines the displayed weights uniquely; its sixth-derivative coefficient is nonzero. The Fourier symbol is $-4h^{-2}{\sin }^2(\theta /2)[1+{\sin }^2(\theta /2)/3]\le 0$. Thus the semidiscrete heat equation $\dot{U}=D_{4,h}U$ is contractive in the discrete L2 norm by the discrete Parseval identity. This says nothing by itself about stability of a subsequently chosen time update.

Take downward displacement $s$ as positive, assume a constant gravitational acceleration $g>0$, and release a point mass from rest at $s=0$. Integrating $\ddot{s}=g$ twice givesFor a downward distance $h>0$, the arrival time and speed areFor the model values $h=20\mathrm{m}$ and $g=9.81\mathrm{m}{\mathrm{s}}^{-2}$, these are approximately $2.0193\mathrm{s}$ and $19.8091\mathrm{m}{\mathrm{s}}^{-1}$. Doubling the height multiplies both quantities by $\sqrt{2}$; it does not double them. The same speed follows from the work-energy theorem, $mgh=m{v}^2/2$.

For intervals of length $\Delta t$, the successive downward displacements areThus equal spacing in time produces increasingly separated positions. This is the odd-number pattern in Galileo's result. It assumes constant acceleration and release from rest; a nonzero initial velocity adds a linear term to $s(t)$.

The approximation can be estimated for a spherical source of radius $R$ and standard gravitational parameter $\mu$. Writing $g_0=\mu /R^2$, the gravitational magnitude at height $h$ isFor $h/R\ll 1$, a Taylor expansion gives $g(R+h)/g_0=1-2h/R+O((h/R{)}^2)$. Small height relative to the source's radius makes nearly constant gravity plausible; negligible drag force remains a separate requirement.

For the stipulated Klein-Gordon packet velocities, take the original normalized spectral amplitudewith $k_0>0$ and $0<\sigma \ll k_0$. Here $\sigma$ is the width of the spectral amplitude, not the standard deviation of its squared magnitude. Initially $\psi (x,0)=e^{-{\sigma }^2{x}^2/2}{e}^{i{k}_{0}x}$, a Gaussian function multiplying a carrier wave. Taking $u=\mathrm{Re}\psi$ gives a real solution of the Klein-Gordon equation with the corresponding initial velocity. The tiny negative-$k$ tail is retained in this exact integral; the statement that it is mainly right-going is a narrow-spectrum approximation, not a claim of strictly positive wavenumbers.

The dimensionless figure uses $c={\omega }_0=k_0=1$ and $\sigma =0.18$ at $t=0,4,8$. Blue curves are $\mathrm{Re}\psi$, orange bounds are $\pm |\psi |$, and the marked carrier wave reference is $x=v_{p}t$; it need not be a crest of the full wavefield. The dashed envelope-centre reference $x=v_{g}t$ is a leading approximation. These are calculated curves, not experimental records. The exact dispersion relation makes the packet slightly asymmetric; a quadratic approximation predicts symmetric broadening about $v_{g}t$.

The first-order Taylor expansion translates the initial wave envelope without changing its shape. Retaining curvature instead givesThe distinction between translation and curvature-driven spreading is the same one exposed by the Gaussian pulse calculation in Pulse Propagation by Richard Fitzpatrick (equations 9.8–9.18). Its first-order approximation is not an assertion that a dispersive packet keeps its shape indefinitely.

Evaluating the quadratic Gaussian integral gives the original approximationand henceHere $L$ is the amplitude width: the standard deviation of the normalized intensity $|{\psi }_2{|}^2$ is $L/\sqrt{2}$. The peak falls as the packet widens, while $\int |{\psi }_2{|}^{2}dx$ stays constant. This conserved integral is not the complete physical wave energy of the Klein-Gordon equation.

Rigid translation requires $|\beta |t{\sigma }^2\ll 1$. The quadratic approximation can allow that quantity to be order one, but still requires the accumulated neglected cubic wave phase to be small over the appreciable spectrum, for example $(t/6){\sup }_{|q|\le 3\sigma }|{\omega }^{\prime \prime \prime }(k_0+q)|(3\sigma {)}^3\ll 1$. This sufficient check controls only the stated three-width spectral window; the Gaussian tails must also be negligible at the required accuracy. Higher derivatives or very long times require returning to the exact spectral integral. The width $L(t)$ is therefore a controlled local approximation; it neither predicts front velocity nor establishes a universal identity between group velocity and information velocity.

Repeated integration by parts gives the Gamma integralThe logarithm of its integrand is $n\log t-t$, whose unique maximum is at $t=n$. Write $t=n+\sqrt{n}s$. For fixed $s$ and large $n$, a Taylor expansion givesThe Big O notation here is for a fixed $s$, or uniformly for $s$ in a fixed bounded interval. It does not claim a uniform approximation over the entire growing integration domain. Consequently the peak, divided by its height, approaches $e^{-s^2/2}$.

The exact substitution givesReplacing the integral by the Gaussian integral ${\int }_{\mathbb{R}}{e}^{-s^2/2}ds=\sqrt{2\pi }$ produces $S_n$. Thus the estimate is peak height $n^n{e}^{-n}$ times width $\sqrt{n}$ times the area of a standard Gaussian peak. This is Laplace's method in a concrete form; controlling the integral's tails justifies the replacement. The next terms, including their products after exponentiation, produce the Stirling correction by Gaussian moments.

For a binary star viewed far from its centre of mass, let $\mathbf{r}$ be the separation and $\mu$ the reduced mass, $M=M_1+M_2$ the total mass, $r=|\mathbf{r}|$ and $x=|\mathbf{x}|$. The Taylor expansion of the point-mass Newtonian gravitational potential has zero dipole because ${\sum }_A{M}_A{\mathbf{x}}_A=0$. Its second moment is ${\sum }_A{M}_A{x}_{Ai}{x}_{Aj}=\mu r_i{r}_j$, soThis holds for an instantaneous separation in any orbit, not only a circular orbit. The mass quadrupole moment convention $Q_{ij}=\mu (r_i{r}_j-r^2{\delta }_{ij}/3)$ gives ${\Phi }_2=-3G{Q}_{ij}{x}_i{x}_j/(2x^5)$.

The group velocity of a differentiable, real dispersion relation branch isIt describes the leading translation of a narrow-band wave packet centred on wavevector ${\mathbf{k}}_0$, assuming a linear model, a locally uniform medium and negligible loss. Narrow-band means that the appreciably excited wavevectors occupy a small neighbourhood on which one smooth branch can be linearized. Write its complex wavefield asThe Taylor expansion $\omega ({\mathbf{k}}_0+\mathbf{q})\simeq {\omega }_0+{\mathbf{c}}_g\cdot \mathbf{q}$ givesThus the rapid wave phase follows the phase velocity, while the slowly varying wave envelope follows the group velocity. This Fourier analysis argument is an outline of the approximation, rather than a statement that every wave packet translates rigidly forever. In one dimension a useful short-time condition is $|{\omega }^{\prime \prime }(k_0)|t(\Delta k{)}^2\ll 1$, where $\Delta k$ is the characteristic spectral width; neglected higher derivatives must also remain small. Curvature of the dispersion relation eventually changes the wave envelope.^[ref]

Two nearby plane waves make the derivative visible without an integral. For equal wave amplitudes $a$, define $\bar{k}=(k_1+k_2)/2$, $\bar{\omega }=({\omega }_1+{\omega }_2)/2$, $\Delta k=k_2-k_1$ and $\Delta \omega ={\omega }_2-{\omega }_1$. The superposition principle gives the exact identityThe periodic wave envelope has speed $\Delta \omega /\Delta k$, approaching group velocity as the two wavenumbers coalesce, while the carrier wave phase has speed $\bar{\omega }/\bar{k}$. Two tones alone form an infinite beat train, not a localized wave packet. The figure shows an exact two-tone surface gravity wave model with dimensionless $\omega =\sqrt{k}$: $k_1=0.9$, $k_2=1.1$. Its carrier phase velocity is approximately $0.99875$, its finite-difference envelope velocity approximately $0.50063$, and ${\omega }^{\prime }(1)=0.5$ is the limiting group velocity. The marked carrier phase line need not coincide exactly with a maximum of the full modulated wave.

The deep-water surface gravity wave therefore has $c_g=c_p/2$, whereas the deep-water capillary wave limit has $c_g=3c_p/2$: crests can drift either forward or backward relative to a wave envelope. For finite-depth surface gravity waves, differentiating the dispersion relation givesHere $\sinh$ is the hyperbolic sine. The shallow-water limit has $c_g\simeq c_p\simeq \sqrt{gh}$; the deep-water limit has $c_g\simeq c_p/2$.^[ref]

In the wave tank, produce a finite train with a smooth turn-on and turn-off and time the same wave envelope feature at two probes separated by distance $L$; its transit estimate is $L/\Delta t$. Independently follow wave crests through that train. Techet's 2003 lab asks students to observe a deep-water train's foremost crest disappearing as the group moves forward, and to compare deep- and shallow-water cases. The interpretation is that a wave crest overtakes the slower wave envelope, dwindling at its front while other crests become visible at its rear. This is a qualitative expected observation in the teaching protocol, not a published set of measured curves.^[ref]

A useful controlled version keeps the paddle wave amplitude small, records $h$ and $\lambda$, and checks both $ka\ll 1$ and $a/h\ll 1$. Compare runs with $kh\gg 1$ and $kh\ll 1$, rather than assuming every ripple tank is shallow relative to every wavelength. Increase train duration at fixed central frequency to narrow its spectrum, and check that the envelope transit estimate stabilizes. Calibrate both probes, sample fast enough to resolve the carrier, and compare records before reflected waves return; an absorbing beach reduces this contamination. A current changes laboratory velocities, large wave amplitudes introduce nonlinear effects, and very short wavelengths require surface tension. Strong damping or noticeable reshaping makes the chosen threshold or peak an unreliable marker of simple group velocity. These are proposed controls and model limits, not claims that the cited lab measured every correction.

Three velocities should consequently be kept separate. The group velocity predicts narrow-band wave envelope motion. An energy-transport velocity is energy flux divided by energy density, using the appropriate time average; it agrees with group velocity in standard lossless linear wave models when all stored energy is included.^[ref] A signal-front velocity concerns the earliest response to a new disturbance and requires a causal model over its full spectrum. A strictly band-limited wave packet has no sharply bounded onset, and a mathematical two-tone envelope cannot test front arrival. In strongly absorbing or broadband propagation, a single group velocity need not describe the observed peak, energy transport, or first detectable response. A complex $\omega$ or $k$ also prevents interpreting its derivative as an ordinary real transport velocity without further analysis.

In particular, an unusually large or negative group velocity from pulse reshaping does not establish superluminal information transmission. For the causal dielectric model analysed by Fitzpatrick, the field ahead of the front vanishes for $t<x/c$, even when its phase velocity exceeds the vacuum speed of light $c$. That result concerns the electromagnetic model's causality; the ideal water-wave diagrams above contain no relativistic signal-front claim.^[ref]

Finally, zero wavenumber requires care: the positive-frequency wave equation branch $c|k|$ is not differentiable at $k=0$, so it has no unique group velocity there, although its one-sided values are $\pm c$. A smooth gapped branch such as the Klein-Gordon equation instead has ${\omega }^{\prime }(0)=0$ despite undefined phase velocity at $k=0$. The same signed branch must be used in both derivatives and ratios; a negative group velocity means an envelope moving in the negative coordinate direction, not a negative elapsed time. The practical value of group velocity is precisely this translation between a local slope on a dispersion diagram and a measurable packet transit, with its assumptions made explicit.

The Gaussian envelope in a Klein-Gordon medium supplies a quantitative next step: the first derivative translates the wave envelope, the second derivative spreads it, and higher derivatives limit that approximation. At finite spectral width, a peak's measured speed need not equal the group velocity at the central wavenumber. Neither a peak nor a fixed fraction of its height generally identifies front velocity.

Here Heun's method means the explicit endpoint-averaging Runge-Kutta method: take an Euler predictor $Y_2=y_n+hf(t_n,y_n)$, sample $k_2=f(t_n+h,Y_2)$ there, and average it with $k_1=f(t_n,y_n)$ for the output. The second Runge-Kutta stage is evaluated at the predictor, not at the corrected endpoint. In particular it requires no implicit solve as the ODE trapezoidal rule does. Naming conventions vary between books; the formulas specify the method unambiguously.^[ref]

For smooth $f$, a Taylor expansion about $(t_n,y_n)$ gives $k_2=f+h(f_t+f_{y}f)+O(h^2)$, where $f_y$ is the state Jacobian matrix. Consequently the output is $y_n+hf+h^2(f_t+f_{y}f)/2+O(h^3)$, matching the exact solution through degree two by the chain rule. This gives local truncation error $O(h^3)$ and, with the usual bounded smoothness and Lipschitz continuity hypotheses on a fixed finite interval, global error $O(h^2)$. The Euler-Heun embedded pair compares this output with its Euler predictor.

For $y^{\prime }=y$, $y_n=1$ and $h=0.2$, the sampled slopes are $1$ and $1.2$. Transporting the second slope to the starting point would end at $1.24$; the average slope $1.1$ instead ends at $1.22$. The exact exponential function reaches $e^{0.2}\approx 1.22140275816$. The magnified endpoint comparison separates the pair difference $0.02$ from the actual Heun local defect $0.00140275816$.

Suppose a function near a critical point has Taylor expansion $f(p+h)=f(p)+P_m(h)+o(\Vert h{\Vert }^m)$, where $P_m$ is its first nonzero homogeneous polynomial term. If $P_m$ takes both positive and negative values on directions, then $p$ is a saddle point of a scalar function. Indeed, on rays $h=rv$ the leading term is $r^m{P}_m(v)$; for each chosen direction its sign dominates the remainder for small positive $r$. This works even when the Hessian matrix is zero. For example, $P_4(\xi ,\eta )=-{\xi }^3(\xi +\eta )$ is negative on $(1,0)$ and positive on $(1,-2)$, hence detects a saddle despite having degree four.

In the Known-variance Gaussian mixture EM example, initialize both means at the same value $m$. The mixture responsibilities are then $1/2$ for every observation. Each M-step sets both means to the sample mean, zero here; the next E-step again gives $1/2$. Thus $(0,0)$ is a fixed point with log-likelihood $-14.593631199$, despite the better fit found from $(-1,1)$.

One can check the failure without trusting an iteration plot. Separate the means as $(-a,a)$ and use the hyperbolic cosine identity:The Taylor expansion uses ${\sum }_i{y}_i^2=18.16$. Small nonzero separation therefore improves the log-likelihood, whereas translating both means to $m$ changes it by $-3m^2$. This stationary point has both improving and worsening directions. Identical initialization cannot discover the separation because every observation keeps identical label probabilities. Trying several distinct initial fits and comparing their final log-likelihoods is useful, but does not prove global optimality.

Keeping the variance fixed also removes a different pathology. With known $\sigma =1$, every component probability density function is at most $(2\pi {)}^{-1/2}$, so $\ell \le -3\log (2\pi )$. If separate component variances are allowed to shrink without restriction, one mixture component can concentrate at an observation while another retains positive density at the other observations, making the likelihood function unbounded. The fixed positive variance in this example prevents that collapse; it does not remove initialization traps.