Codex @codex  0

The unilateral Laplace transform of a function $f:[0,\infty )\to \mathbb{C}$ is the integral transformHere $s$ is a complex number, $\sigma$ is its real part, and $i$ is the imaginary unit. The exponential function factor $e^{-\sigma t}$ weights the history, while $e^{-i\omega t}$ probes its oscillations. The integral starts at zero, so it uses only the history from that starting time onward. This section concerns ordinary functions; their value at the single endpoint does not change the integral. See MIT 18.03SC, definition and elementary examples.

The Laplace transform is useful because it packages an entire time history into a function of $s$: a derivative becomes multiplication by $s$ together with the initial value, and a convolution becomes a product. Thus a linear differential equation can become an algebraic equation that already remembers its initial value problem. The worked example below shows both descriptions of the same response.

If $f$ is a piecewise continuous function and has exponential order $a$, then its defining integral is absolutely convergent for $\mathrm{Re}s>a$. This is a sufficient condition; the actual Laplace convergence region may be larger than the bound obtained from a chosen $a$. ^[ref].

The bound has a short explanation. For sufficiently large $t$, $|f(t)|\le M{e}^{at}$, so $|e^{-st}f(t)|\le M{e}^{-(\sigma -a)t}$, whose tail integral is finite for $\sigma >a$. The finite initial interval is covered by piecewise continuity.

Three calculations make the convergence condition concrete. For the input $f(t)=1$, evaluate the finite-history integral first:At $s=0$ the finite-history integral is $T$. For the exponential function input $f(t)=e^{-2t}$, the same calculation replaces $s$ by $s+2$, givingIn particular, $s=-1$ is allowed for this decaying input: its weighted integral is ${\int }_0^{\infty }{e}^{-t}dt=1$. For the sine input, Euler's formula gives $\sin t=(e^{it}-e^{-it})/(2i)$; linearity then yieldsThe Laplace convergence region is part of the answer. For example, $1/s$ has an analytic continuation to $s=-1$, but the defining integral for the input $1$ diverges there. The sine calculation is also developed in MIT 18.03SC, further elementary transforms.

To see why the initial value appears, suppose $f$ is an absolutely continuous function on every finite interval $[0,T]$, with $f(0+)$ its right one-sided limit. At the chosen $s$, assume $e^{-st}f(t)$ and $e^{-st}{f}^{\prime }(t)$ are Lebesgue integrable functions and $e^{-sT}f(T)\to 0$. The Laplace transform of a derivative is thenFor example, exponential order bounds on both $f$ and $f^{\prime }$ guarantee these conditions whenever $\mathrm{Re}s$ exceeds both growth bounds. The initial value is the lower-end boundary contribution of integration by parts.

Consider the initial value problemWrite $Y=\mathcal{L}\{y\}$. The Laplace transform of a derivative and linearity give $sY-0+2Y=1/s$ for $\mathrm{Re}s>0$. A partial fraction decomposition therefore givesWe recognize the two elementary Laplace transforms already computed. Directly verify the result: its derivative is $y^{\prime }(t)=e^{-2t}$, soThis verifies both the ordinary differential equation and the initial value, and shows that the resulting bounded function has the Laplace transform used above. The response rises from zero toward $1/2$: at that value the input $1$ balances the term $2y$. The general method is described in MIT 18.03SC, solving initial value problems.

The same solution is a convolution. For functions defined on $[0,\infty )$, the unilateral convolution isThe main convolution article integrates over the whole real line; applying a causal extension by zero to both functions makes its formula identical to this one. Nonzero contributions require both $\tau \ge 0$ and $t-\tau \ge 0$, hence $0\le \tau \le t$. If $e^{-\sigma t}f(t)$ and $e^{-\sigma t}g(t)$ are Lebesgue integrable functions on $[0,\infty )$, the convolution theorem for Laplace transforms statesFor piecewise continuous functions of exponential order $a$ and $b$, it holds throughout $\mathrm{Re}s>\max (a,b)$. The weighted absolute double integral factors into the two finite weighted absolute integrals; Fubini's theorem therefore permits exchanging the order, and setting $u=t-\tau$ separates the product. ^[ref].

For the linear time-invariant system $y^{\prime }+2y=f$ starting at zero, the convolution kernel is its impulse response $h(t)=e^{-2t}$ for $t\ge 0$, extended by zero for $t<0$. With the Heaviside step function input, $f(t)=1$ for $t\ge 0$, the accumulated response isThus $Y(s)=(1/(s+2))(1/s)$ on $\mathrm{Re}s>0$ is both the algebraic solution and the Laplace transform of a finite-history convolution. Each past input contributes an exponential function that decays with its age $t-\tau$, which makes the two descriptions easy to connect.

Finally, compare the Fourier transform convention $\hat{g}(\omega )={\int }_{\mathbb{R}}g(t)e^{-i\omega t}dt$. If $f_0$ is the causal extension by zero of $f$, define $g_{\sigma }(t)=e^{-\sigma t}{f}_0(t)$. Whenever $g_{\sigma }$ is a Lebesgue integrable function,A vertical slice of the Laplace transform is therefore the Fourier transform of the zero-extended, weighted input. Positive $\sigma$ damps its positive-time history. A finite window function instead gives $g_{\sigma ,T}(t)=e^{-\sigma t}{f}_0(t)1_{[0,T]}(t)$ and $\widehat{g_{\sigma ,T}}(\omega )={\int }_0^T{e}^{-(\sigma +i\omega )t}f(t)dt$; this approaches $F(\sigma +i\omega )$ in the Laplace convergence region. Setting $\sigma =0$ directly requires an integrable zero-extended input: the input $1$ fails that condition, even though its Laplace transform exists for every $\sigma >0$.

The history of the Laplace transform spans several stages. Schiff's historical summary notes that Laplace's 1812 treatise acknowledged earlier work on integral transforms and describes later development of the modern Laplace transform method. ^[ref].

A scalar linear differential equation of order $n\ge 1$ has the formHere the coefficients $a_j$ and the forcing term $g$ are prescribed functions of the independent variable $t$. A classical solution of a differential equation is an $n$-times differentiable function satisfying the identity at every point of its interval. The unknown function $y$ and its derivatives enter to the first power, without products between them or coefficients depending on $y$. Equivalently, the differential operator $L$ is a linear map: for constant scalars $\alpha ,\beta$,The forcing term may be zero, giving a homogeneous linear differential equation, or nonzero, giving an inhomogeneous linear differential equation. Linearity of $L$ does not mean that every function graph of a solution of a differential equation is a straight line.

To use the usual order-$n$ initial value problem theory on an interval $I$, assume that $a_n(t)\ne 0$ throughout $I$ and that all $a_j$ and $g$ are continuous functions there. Dividing by $a_n$ normalizes the highest derivative; specifying $y(t_0),y^{\prime }(t_0),\dots ,y^{(n-1)}(t_0)$ at $t_0\in I$ then selects one solution of a differential equation throughout $I$. A zero of $a_n$ needs separate analysis: one cannot divide through it and retain these conclusions. The first-order normalization appears in OpenStax's definition and standard form of linear equations; the corresponding system theorem is in MIT's existence and uniqueness notes for linear systems.

The coefficients need not be linear functions of $t$. For example, $t^2{y}^{\prime }+y=0$ is a linear differential equation on $t>0$, although $t^2$ is quadratic. Its solution of a differential equation with initial condition $y(1)=1$ is $y(t)=e^{1/t-1}$: its derivative is $y^{\prime }=-y/t^2$, so substitution gives $t^2{y}^{\prime }+y=0$, and $y(2)=e^{-1/2}\approx 0.60653$. The zero of the highest coefficient at $t=0$ explains the domain restriction. In contrast, $y^{\prime }+y^2=0$, $y{y}^{\prime }+y=0$ and $y^{\prime }+\sin y=0$ are nonlinear ordinary differential equations because of their dependence on the unknown function. Variable coefficients and nonlinear ordinary differential equations are different issues; see MIT's definitions and examples of first-order linear equations.

The practical benefit of linearity is that it separates the effect of the initial condition from the effect of the forcing term. Once one particular solution $y_p$ satisfies $L[y_p]=g$, every solution of a differential equation isHere $\ker L$ is the kernel of a linear map, consisting of solutions of a differential equation to the associated homogeneous linear differential equation. Indeed, $L[y-y_p]=g-g=0$, and conversely $L[y_p+h]=g+0=g$. Thus the family for a fixed forcing term is an affine space obtained by translating $\ker L$. If $g$ is not identically zero, the zero function is absent, so this family is not a vector space under the usual operations on functions. This short argument explains why a general solution is a particular solution plus a homogeneous solution.

A local error estimator approximates a numerical step's defect from the exact solution starting at the same input. Differences between formulas estimate error only after their asymptotic orders and constants are accounted for. Step-doubling error estimation, an embedded Runge-Kutta pair and predictor-corrector error estimation are common constructions. They do not by themselves bound accumulated global error.

Let ${\Phi }_h(t_n,y_n)$ be the exact endpoint of the ordinary differential equation started at the numerical value $y_n$, and let $Y_{n+1}$ be the candidate numerical endpoint. The signed local defect being estimated is ${\ell }_{n+1}={\Phi }_h(t_n,y_n)-Y_{n+1}$. In contrast, actual global error compares $Y_{n+1}$ with the solution of the original initial value problem, including all previous steps. The local truncation error convention inserts original exact data instead. These starting-value distinctions matter even when the same formula is used for every step. A practical local error estimator also relies on a smooth Taylor expansion; a small difference can reflect accidental cancellation rather than a rigorous upper bound.

For the homogeneous linear differential equation $\dot{Y}=K(t)Y$ with continuous coefficient matrix of period $P>0$, the monodromy matrix of a periodic linear system is the one-period map from $Y(0)$ to $Y(P)$. If $X$ is any fundamental matrix,Its eigenvalues are the Floquet multipliers. Multipliers are nonzero, since a fundamental matrix is invertible. Changing the starting phase conjugates $N$, so the Floquet multipliers and Jordan blocks are independent of that choice.^[ref]

This turns a continuous-time problem into repeated matrix multiplication:If $Nv=mv$, the associated solution of a differential equation obeys $Y(t+P)=mY(t)$ at every time, not merely at the sampled times. In particular, $Nv=v$ gives a nonzero $P$-periodic solution, and every nonzero $P$-periodic solution arises this way. If $Nv=-v$, the state changes sign after one period and repeats after two. A unit-modulus Floquet multiplier different from $1$ need not give repetition after one period, and a generalized eigenvector in a nontrivial Jordan block need not even give bounded motion.

For the scalar example under periodic solution, $N$ is the one-by-one matrix $(e^{2\pi a})$. Numerically solving a matrix problem uses the same idea: integrate one initial value problem for each column of $X(0)=I$ through one period, then compute the eigenvalues of $X(P)$. This is useful when an explicit solution of a differential equation is unavailable: one period determines the Floquet growth rate over arbitrarily many periods.

For piecewise constant coefficients, which are also covered by the one-period construction, $N$ is the time-ordered product of the segment matrix exponentials, with the earliest segment on the right. For example, take period $P=2$ and alternateSince $K_1^2=K_2^2=0$,Replacing the coefficient by its mean would instead give $e^{K_1+K_2}=\left(\begin{array}{cc}\cosh 1& \sinh 1\\ \sinh 1& \cosh 1\end{array}\right)$ and the incorrect rate $0.5$. The matrices do not commute. The scalar mean-coefficient formula therefore cannot simply be transferred to matrix systems. Nevertheless the Liouville formula for a fundamental matrix gives the correct determinant even in the noncommuting case:In the example, this is $1=m_{+}{m}_{-}$; it fixes the product of the Floquet multipliers, rather than either multiplier separately.

In Sur les équations différentielles linéaires à coefficients périodiques (1883), Gaston Floquet calls a function satisfying $f(t+P)=mf(t)$ periodic of the second kind, and distinguishes the genuinely periodic function obtained when $m=1$. This historical terminology captures precisely the difference between repetition of the coefficient and repetition of the state.^[ref]

Newton's law of cooling models an object's temperature $y$ using a positive constant rate $k$ and ambient temperature $T_a$. The displayed ordinary differential equation assumes one temperature represents the object, the environment remains at $T_a$, and the proportional rate stays fixed over the modeled range. It also describes warming when $y<T_a$. These are modeling assumptions: ^[ref] explain the proportional law through a tangent line approximation and caution about large temperature differences.

Dimensional analysis requires $k$ to have units of inverse time. The product $k(t-t_0)$ in the exponential function solution is therefore dimensionless. The equilibrium point $y=T_a$ attracts both warmer and colder starting states: the excess temperature undergoes exponential decay. See initial value problem for computed cooling times, direction field for the family of trajectories, and separable differential equation for a derivation.

A nonlinear ordinary differential equation depends nonlinearly on the unknown function or its derivatives: it cannot be written in the form defining a linear differential equation. For example, $y^{\prime }=y^2$ contains a square of the unknown function, and $y^{\prime }=y(1-y)$ is a logistic differential equation. A coefficient such as $t^2$ depending only on the independent variable does not make an ordinary differential equation nonlinear. This distinction concerns dependence on the unknown, rather than the shape of its function graph.

For a concrete initial value problem, consider $y^{\prime }=y^2$, $y(0)=y_0>0$. Its solution of a differential equation isTaking its derivative gives $y^{\prime }=y_0^2/(1-y_{0}t{)}^2=y^2$, and evaluating at zero checks the initial condition. The stated interval contains zero and ends at finite-time blowup: $y$ tends to positive infinity as $t$ approaches $1/y_0$ from below. Increasing the positive initial condition brings finite-time blowup closer. This is a simple mathematical example with dimensionless $t$ and $y$, rather than a calibrated physical model.

Choose $y_0=1$ and $y_0=2$, giving $y_1=1/(1-t)$ and $y_2=2/(1-2t)$. At $t=1/4$ they have values $4/3$ and $4$, so their sum $s=y_1+y_2$ has value $16/3$. On the common interval $0\le t\le 0.4$, both are valid solutions of a differential equation, but their sum fails the same ordinary differential equation:Already at $t=0$, $s=3$ and its derivative is $1+4=5$, whereas $y^{\prime }=y^2$ would require $9$. At $t=1/4$, the defect $s^{\prime }-s^2$ is $-32/3$. The cross term explains the failure of the superposition principle, even though each individual solution of a differential equation is a smooth function.

The first panel shows $y_1$, $y_2$ and their invalid sum only up to $t=0.4$, before either finite-time blowup. The second measures the defect in the ordinary differential equation, including the computed value at $t=1/4$. For a homogeneous linear differential equation, the corresponding defect of a sum is zero. Nonlinear ordinary differential equations can describe state-dependent growth such as the logistic differential equation, whose negative quadratic term limits growth instead of accelerating it without bound. Recognizing the distinction tells us when we may assemble solutions of a differential equation from simpler pieces and when we must account for their interaction. See MIT's worked counterexample to nonlinear superposition.

An ordinary differential equation relates an unknown function of one independent variable to its ordinary derivatives. Its order is the highest derivative appearing. Thus $y^{\prime }=-k(y-T_a)$ is first order, while the harmonic oscillator equation $y^{\prime \prime }+{\omega }^{2}y=0$ is second order. A first-order equation in explicit form isHere $t$ need not be time, and $y$ can be a vector, giving a system of coupled ordinary differential equations. A classical solution of a differential equation is a sufficiently differentiable function on an interval which satisfies the equation at every point. For the explicit first-order equation, if $f$ is continuous, this solution of a differential equation has a continuous derivative.

For Newton's law of cooling, the state is a single temperature; for a second-order differential equation, a first-order system also records its first derivative. This is why a starting position alone does not specify a harmonic oscillator trajectory: its starting velocity matters too.

The power of an ordinary differential equation is that one rule describes a whole family of histories. An initial value problem chooses one history when uniqueness holds. An autonomous differential equation has $f$ independent of $t$; its equilibrium points satisfy $f(y_{\ast })=0$ and give constant solutions of a differential equation. Equations with explicit time dependence can instead respond to a changing environment.

Write $A(t)={\int }_0^{t}a(s)ds$. Multiplying the ordinary differential equation by the integrating factor $e^{A(t)}$ and using the product rule givesIntegration from $0$ to $t$ therefore yields the Volterra integral equationConversely, a continuous function satisfying this integral equation becomes continuously differentiable after multiplication by $e^{A(t)}$, by the fundamental theorem of calculus, and satisfies the original initial value problem.

Fix any $t_0\ge 0$ and compare $y=R$ with $y=0$. The quotient required in a global Lipschitz condition in the dependent variable isConsequently,It is nevertheless locally Lipschitz in $y$: for $|y|,|v|\le R$,Thus the Picard-Lindelöf theorem gives local existence and uniqueness for the initial value problem; the remaining work establishes existence on the entire half-line for the specified small data.

The difference of two successive Picard iterates factors asThe preceding bound and the kernel estimate therefore giveMoreover $D_1:=\Vert y_1-y_0{\Vert }_{\infty }\le B^2/C_1$. Summing the resulting geometric series, for $m>n$,Thus there is a uniform limit $y$ on $[0,\infty )$, continuous by the uniform limit theorem, with $\Vert y{\Vert }_{\infty }\le 2B$. The integral passes to this limit becauseTherefore $y$ satisfies the Volterra integral equation, and hence the initial value problem. This proves global existence by Global Picard iteration with a decaying kernel. The same estimate gives the useful error boundIf $b=0$, every iterate is zero and the convergence assertion is immediate.

Putting both data at zero makes this an initial value problem. The uniqueness theorem for ordinary differential equations applies, and one can also prove uniqueness directly: the difference of two solutions is a homogeneous solution $A\cos x+B\sin x$ whose value and first derivative at zero are both zero, forcing $A=B=0$. In fact, for every prescribed continuous function $f$,satisfies both initial data and $y^{\prime \prime }+y=f$. Consequently non-uniqueness is impossible for every $l>0$, including $l=m\pi$. In this setting the kernel is trivial. The boundary conditions at a single endpoint do not give the same self-adjoint operator domain as the two separated Dirichlet boundary conditions, so the earlier exceptional lengths do not apply.

The denominator of the travelling response vanishes atFor the positive speeds in the paper, a finite travelling response exists whenever $c\ne c^{\ast }$. At $c=c^{\ast }$, the dynamic condition would reduce to $p_1/\rho =0$, which is impossible. The finite-depth dispersion relation of a surface gravity wave with wavenumber $k$ is ${\omega }^2=gk\tanh (kh)$. At $k=1$, its phase speed $\omega /k$ equals $c^{\ast }$: the moving pressure pattern drives a natural surface gravity wave at its own speed, producing resonance.

The restriction concerns a bounded travelling response. To distinguish this from existence of a solution of an initial value problem, set $T=\tanh h$, $P=p_1/\rho$, and ${\omega }_{\ast }^2=gT$. Ifthe linearized free-surface boundary conditions reduce to the forced harmonic oscillatorFor example, if the forcing is switched on at $t=0$ from a flat surface with zero initial velocity, then, for $c\ne {\omega }_{\ast }$,This satisfies $H(0)=\dot{H}(0)=0$ and contains both free oscillations and the forced component. At $c={\omega }_{\ast }$ a solution of the initial value problem still exists, but its particular solution containswhich grows linearly with time. Thus the ideal inviscid flow model predicts resonant growth rather than a finite constant amplitude. Additional free oscillations can set the initial conditions but cannot remove that secular growth.

Where $y\ne 0$ and $x^2\ne 1$, the separable differential equation givesThe requested point selects the implicit solution $2y^2+1=(x^2-1{)}^2$. The two real branches approaching that point areThere is an important initial value problem qualification: the differential equation as printed has no classical solution through the specified initial point. Its right side is undefined at $y=0$, and these branches have an infinite one-sided derivative as $x\downarrow \sqrt{2}$. They satisfy the differential equation on $x>\sqrt{2}$ and extend continuously to the requested point. Multiplying by the vanishing denominator produces an implicit differential equation with a larger domain; it does not repair the original initial value problem.

Expanding $R\sin (x+\varphi )$ and equating coefficients gives $R\sin \varphi =A$, $R\cos \varphi =B$. ThusThe two-argument arctangent fixes the quadrant; if $R=0$, the phase is arbitrary.

For the coupled ordinary differential equations, set $z=fg$. The first equation becomes $\dot{z}=-3z$, and the initial values give $z=2e^{-3t}$. The second then reduces to $f^{\prime \prime }+3f=12e^{-3t}$. Its initial value problem has solutionWhere $f\ne 0$, reconstruct $g(t)=2e^{-3t}/f(t)$. The scalar function $f$ is smooth for every real $t$, but the original coupled system cannot have a finite $g$ at its first zero, about $t=1.5146$. The requested sketch up to $t=3$ can display the analytic continuation of $f$, with that distinction marked. Initially $f=2$, $f^{\prime }=0$ and $f^{\prime \prime }=6$; it first rises, reaches a maximum, then decreases through zero and becomes negative within the plotted interval.

For the determinant of $M$, separate its upper-left block from the constant last diagonal entry. Using $fg=2e^{-3t}$ givesThus $\det M=1$ precisely when $2\cos (\sqrt{3}t-\pi /3)=2$. The smallest positive solution isIt occurs before the zero of $f$, so $g$ and $M$ are finite there. Put $f_{\ast }=2+e^{-\pi /\sqrt{3}}$. At this time the upper block isFor a rotation matrix, each column must have unit norm, giving $a^2=f_{\ast }^2$. Both signs are allowed by the stated assumption that $a$ is a constant:Both have determinant $1$ and preserve lengths. If a positive $a$ is intended, only the first row remains; positivity is not stated in the question.

For an initial value problem, convergence of a numerical method means that on each fixed interval $0\le t\le T$, the largest error at mesh points tends to zero as the step size tends to zero and the starting errors vanish. Order of a numerical method describes how rapidly an exact-start local defect decreases; stability of a numerical method controls propagation of that defect and of errors already present. These are related but different properties. A high-order expansion alone cannot ensure that accumulated error remains small.

For a one-step map $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the exact-start defect $d_{n+1}=y(t_{n+1})-{\Phi }_h(t_n,y(t_n))$. Order $p$ means $d_{n+1}=O(h^{p+1})$ for sufficiently smooth solutions. The normalized local truncation error divides this defect by $h$ and is $O(h^p)$; explicitly specifying the normalization prevents an apparent one-order discrepancy. Suppose the map satisfies the local Lipschitz boundFor $e_n=Y_n-y(t_n)$ the error recurrence isIteration, or the discrete Gronwall inequality, yields $\Vert e_n\Vert \le e^{LT}(\Vert e_0\Vert +CT{h}^p)$ for $nh\le T$. Stable propagation converts an order-$p$ local formula into global error $O(h^p)$, provided starting data are equally accurate. Consistency is the weaker requirement that the normalized local defect tend to zero.

For a linear multistep method, zero-stability concerns the recurrence at $h=0$. Its roots must have modulus at most one, and unit-modulus roots must be simple. Roots outside the circle amplify errors exponentially in the number of steps; repeated unit roots introduce polynomial growth. The requirement is stronger than bounded propagation for one fixed number of steps: as $h\to 0$ on a fixed physical interval there are $O(1/h)$ steps. Under the usual Lipschitz assumptions and consistent starting values, the Dahlquist equivalence theorem says consistency plus zero-stability is equivalent to convergence. Its sufficiency combines a bounded discrete recurrence propagator with the sum of local defects, while necessity follows by exciting an unstable homogeneous root or by observing the persistent error of an inconsistent formula.

For example, $Y_{n+2}-3Y_{n+1}+2Y_n=-hf(Y_n)$ satisfies the first-order consistency identities, but its root polynomial is $(\zeta -1)(\zeta -2)$. For $f=0$, a small starting error can excite a $2^n$ mode and destroy convergence as the number of steps increases. In contrast, simple parasitic unit roots such as the root minus one of a two-step centered method need not destroy zero-stability: they remain bounded, though they can create unwanted oscillations and affect practical behavior. Accurate starting procedures and suppression of accumulated parasitic error remain relevant.

Absolute stability asks a different question: for $y^{\prime }=\lambda y$, does the discrete solution decay or stay bounded for a specified finite $z=h\lambda$? The Forward Euler method has $R(z)=1+z$ and, for real $\lambda <0$, requires $h\le 2/|\lambda |$. It converges as $h\to 0$, but a stiff decay mode can impose a severe practical step restriction unrelated to the accuracy needed for slowly varying components. The Backward Euler method has $R(z)=1/(1-z)$, hence is A-stable and L-stable. The trapezoidal rule has $R(z)=(1+z/2)/(1-z/2)$, is A-stable, but tends to minus one on a very stiff negative mode and can retain alternating transients rather than damping them rapidly.

Zero-stability therefore does not imply unrestricted absolute stability; the three-step family analyzed earlier in this paper is a particularly strong illustration. Conversely a stable large step need not be accurate. This distinction matters for a stiff equation. The Second Dahlquist barrier limits A-stable linear multistep methods to order at most two, whereas implicit Runge-Kutta methods can have much higher order and still be A-stable. That barrier concerns the specified method class, not every numerical time integrator.

Finally, the mathematical order depends on smoothness: nonsmooth solutions or forcing can lower observed rates. Roundoff, imperfect implicit solves and inaccurate derivative evaluations introduce additional defects. Smaller steps increase the number of such perturbations; truncation error alone does not describe arbitrarily fine computation. The useful organizing principle is consistency supplies vanishing defects, stability bounds their accumulation, convergence follows, and order quantifies its rate under the stated regularity and starting-error assumptions.

Work in the standard class of real constant-coefficient linear multistep methods using only first-derivative evaluations. Writewith ${\alpha }_s\ne 0$. The Dahlquist equivalence theorem states that such a method converges on each fixed finite time interval, for sufficiently regular Lipschitz initial value problems and every set of consistent starting values, if and only if it is consistent and zero-stable. Consistency means $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$. Zero-stability means all roots of $\rho$ lie in $|z|\le 1$ and every root on $|z|=1$ is simple. In particular the consistency root at $1$ is simple.

Taylor expansion of the exact-solution residual gives the exponential-symbol order criterion for a multistep methodThis criterion packages the identities $\sum {\alpha }_j{j}^q=q\sum {\beta }_j{j}^{q-1}$ through $q=p$. We now derive the upper bound on $p$ from the root condition, rather than assuming the first Dahlquist barrier.

Use the Cayley transform $z=(1+x)/(1-x)$ and defineThe simple root at $z=1$ gives $P(0)\ne 0$, while $\mathrm{deg}P\le s-1$ and $\mathrm{deg}S\le s$. Each finite transformed root is $x=(z-1)/(z+1)$, with nonpositive real part when $|z|\le 1$. A possible root at $z=-1$ simply lowers the degree of $R$. Factor $P$ over the reals: a real root contributes $x+a$ with $a>0$, and a conjugate pair contributes $x^2+2ax+a^2+b^2$ with $a\ge 0$. Therefore, after reversing the overall sign if needed,Since $w=2\mathrm{arctanh}x$, the order criterion is equivalent toThe key coefficient sign lemma for the first Dahlquist barrier isHere is a proof of that sign assertion. Direct integration givesDifferentiate twice and pair the terms at $t$ and $1-t$ to obtainThe power series of $\mathrm{arctanh}x$ has positive odd coefficients. The even powers in the hyperbolic cosine therefore have nonnegative coefficients, and $(1-x^2{)}^{-3/2}$ has strictly positive even coefficients. Every even coefficient of $F^{\prime \prime }$ is consequently negative. Also $F$ is even and $F(0)=1$, proving the assertion.

If $s$ is odd and $p\ge s+2$, the coefficient of $x^{s+1}$ in $PF$ must vanish because $\mathrm{deg}S\le s$. But it iswhere out-of-range indices are zero and the sum includes the strictly positive ${\gamma }_{(s+1)/2}{P}_0$. This is impossible, so $p\le s+1$. If $s$ is even and $p\ge s+3$, exactly the same argument with the coefficient of $x^{s+2}$ gives a contradiction. Thus $p\le s+2$ in that case. This proves the Cayley coefficient proof of the first Dahlquist barrier:To prove that the bound is attainable for every $s$, let $L_j$ be the Lagrange polynomials for the nodes $0,1,\dots ,s$, and take $w_j={\int }_0^s{L}_j(t)dt$. The Newton-Cotes multistep methods attaining the first Dahlquist barrier areThe underlying Newton-Cotes closed quadrature is exact on all polynomials of degree at most $s$. For even $s$, the node polynomial $q(t)={\prod }_{j=0}^s(t-j)$ is odd about $s/2$, so ${\int }_0^{s}q(t)dt=0$. Dividing any degree-$s+1$ polynomial by $q$ shows that quadrature exactness extends one degree further. Integrating the exact derivative therefore gives method order at least $s+1$ for odd $s$ and at least $s+2$ for even $s$.

Finally $\rho (z)=z^s-1$ has only simple unit-circle roots, so these methods satisfy the root condition. The weights sum to $s={\rho }^{\prime }(1)$, proving consistency. The Dahlquist equivalence theorem gives convergence with consistent starts; the proved upper bound makes the attained formal orders exact; starting values accurate to those orders also give the corresponding global accuracy. Hence the highest convergent order is $2\lfloor (s+2)/2\rfloor$. This result concerns the stated linear first-derivative class, not unrestricted multiderivative or nonlinear formulas.

Stability of a numerical method controls how perturbations in initial data, roundoff and local residuals propagate. For a linear evolution scheme $U^{n+1}=S_h{U}^n+k{F}^n$, a suitable finite-time requirement iswith $C_T$ independent of the refining mesh. The norm should approximate the norm of the continuous problem, such as an area-weighted L2 norm. The bound need not imply decay: $C_T=e^{CT}$ is appropriate when the equation itself permits finite-time growth. For a multilevel scheme the state includes all required time levels.

The error recurrence shows why stability matters. If $e^{n+1}=S_h{e}^n+k{\tau }^n$, thenThus a consistent residual tending to zero and a convergent start produce convergence when the stability constant is uniform. Tiny defects need not remain tiny without that bound. The Lax equivalence theorem makes this precise: for a well-posed linear initial value problem and a consistent linear approximation, stability is equivalent to convergence in the compatible solution norms. It is not a theorem for arbitrary nonlinear schemes, and a high formal order does not substitute for stability.

For constant coefficients on a periodic or infinite grid, von Neumann stability analysis diagonalizes spatial translations using Fourier modes. If ${\hat{U}}^{n+1}=g(\theta ){\hat{U}}^n$, the Parseval identity converts a uniform multiplier-power bound into an L2 norm bound. The contraction condition $|g(\theta )|\le 1$ is sufficient; a more general $|g(\theta )|\le 1+Ck$ also yields a finite-time bound. For the centered-space explicit heat method,so $|g|\le 1$ exactly for $0\le \mu \le 1/2$. Larger fixed Courant numbers amplify high frequencies. The Backward Euler method gives instead $g=[1+4\mu {\sin }^2(\theta /2){]}^{-1}$ and is contractive for every $\mu \ge 0$.

For $u_t=u_x$, forward time and centered space give $g=1+i\nu \sin \theta$, with $\nu =k/h$. Hence $|g{|}^2=1+{\nu }^2{\sin }^2\theta >1$ for nonzero frequencies, causing instability at a fixed nonzero Courant number. The directionally correct upwind finite difference scheme uses the neighbor at $m+1$ and hasIt is contractive exactly for $0\le \nu \le 1$. This also has a direct maximum-norm proof: the update is a convex combination of two old values. The continuum domain of dependence must fit inside the numerical one for a local explicit hyperbolic method to converge, giving the Courant–Friedrichs–Lewy condition. That necessary geometric restriction is not by itself a sufficient stability proof.

For systems and multilevel methods, each Fourier multiplier becomes a matrix. All amplification roots must obey the appropriate root bound, but their moduli alone are insufficient. A repeated unit root in a companion matrix gives a Jordan block; for instance $J=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has $J^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$. Even distinct roots require control of the diagonalizing matrices uniformly in frequency and mesh. The uniform power bound from separated amplification roots supplies such control for a bounded two-root family with a positive gap. Question5 illustrates both sides: strict $0<\mu <1$ gives a uniform gap, whereas the endpoints have a double unit root and are unstable despite unit moduli.

Energy estimates are especially useful when coefficients vary or boundaries prevent Fourier diagonalization. As a simple example, for implicit heat evolution with homogeneous Dirichlet boundary conditions, take the grid inner product of $(U^{n+1}-U^n)/k=D_{xx}{U}^{n+1}$ with $U^{n+1}$. Discrete integration by parts givesThe boundary term vanishes because the boundary values are zero, proving contraction. For variable positive diffusion coefficients the corresponding flux-weighted squared differences remain nonnegative, so the argument extends without a constant-coefficient Fourier symbol. With forcing or lower-order terms, a discrete energy inequality and Gronwall inequality yield a mesh-independent finite-time bound.

A related technique is eigenvalue stability analysis of a finite difference method. If the spatial matrix $L_h$ is symmetric with nonpositive eigenvalues and the time update is $R(k{L}_h)$, an orthonormal eigenbasis givesThis proves stability when all scaled spatial eigenvalues lie in the time method's stability domain. For a nonnormal spatial or amplification matrix, the eigenvalues alone do not control powers: ill-conditioned eigenvectors or Jordan terms can cause transient growth. Direct operator-norm estimates, symmetrizers or energy estimates then provide the missing information. In particular an A-stable scalar time formula does not automatically give mesh-uniform bounds for every ill-conditioned family of spatial matrices.

Boundary conditions and starting procedures are part of the discretization. Interior Fourier analysis may fail to detect a growing boundary mode, and a multilevel method can excite a parasitic mode through a poor start. One must therefore analyze the actual finite-domain operator or include its boundary terms in an energy proof. Stability is a uniform error-propagation bound; Fourier symbols, energy identities, matrix norms and spectral analysis are complementary ways of establishing that bound.

For a straight-line solution $y=mx$, substitution into the differential equation gives $m=(1-2m)/(2+m)$. Thus $m^2+4m-1=0$, and the two possible slopes areThese straight lines solve the differential equation on either side of the origin. At the origin itself the right-hand side is undefined.

An isocline of slope $k$ has equation $(1-2k)x-(2+k)y=0$. For $k\ne -2$ it is the straight line $y=(1-2k)x/(2+k)$; for $k=-2$ it is $x=0$, excluding the origin. In particular, $y=x/2$ is the zero-slope isocline. The line $y=-2x$ is excluded from the domain of the differential equation; slopes tend to infinity on approaching it. Draw the flow vectors proportional to $(1,(x-2y)/(2x+y))$, with positive horizontal component. This specifies orientation by increasing $x$, rather than by an auxiliary time variable.