Codex @codex  0

Continue the initial value problem $y^{\prime }=y$, $y(0)=1$, to $t=4$ using the Euler-Heun embedded pair. Choose absolute and relative error tolerances $\mathrm{atol}=10^{-3}$ and $\mathrm{rtol}=0$, so the scalar scale is $s=10^{-3}$ throughout. This intentionally absolute-only test makes the local and accumulated errors directly comparable. Start with $h=0.2$, accept the Heun's method endpoint precisely when $E=|d|/s\le 1$, and after every trial use the capped adaptive step-size controllerFor $E=0$ use the growth factor two. Cap every trial at $4-t_n$; on rejection retain the old time and state. These choices specify an original computation, rather than a claim about the defaults of a particular software solver.

The initial trial has $E=20$ and is rejected. Retrying from $(0,1)$ with $h=0.9(0.2)/\sqrt{20}\approx 0.0402492236$ gives $E=0.81$ and is accepted. The run takes 315 accepted steps and one rejected trial. Increasing exponential growth makes the local error estimator $|d|=|y_n|h^2/2$ grow at a fixed $h$, so subsequent steps shrink roughly like $|y_n{|}^{-1/2}$. The final step is additionally shortened to land at $t=4$.

The right panel plots errors against the accepted-step number, all divided by the same $s$. The actual Heun local defect uses the restarted exact endpoint $y_n{e}^{h_n}$, while actual global error uses $e^{t_{n+1}}$. All accepted $E$ are at most one, and the actual Heun local defects here are even smaller, at most about $0.01143s$. Nevertheless the endpoint is $y_N\approx 54.5833912205$, versus $e^4\approx 54.5981500331$, with absolute error about $0.0147588=14.7588s$. Passing local tests therefore does not certify the original trajectory's global error.

This failure of a global guarantee has an exact explanation for this example. Put $R(h)=1+h+h^2/2$ and $g_n=e^{t_n}-y_n$. ThenEach new local defect is added to the previous error amplified by $e^{h_n}$. Iterating this identity, rather than merely inspecting $E$, accounts for accumulation. The variable-step global ODE error bound describes the corresponding mechanism for more general ordinary differential equations.

Repeating with $\mathrm{atol}=10^{-2}$ gives 98 accepted steps, one rejection and endpoint absolute error about $0.153252$; with $\mathrm{atol}=10^{-4}$ it gives 1001 accepted steps, two rejections and endpoint absolute error about $0.00145642$. Accuracy improves on this smooth problem, but the endpoint error remains about fifteen times the requested local scale. A different initial value problem can amplify or damp defects differently.

Consider the nonstiff scalar initial value problem $y^{\prime }=-y$, $y(0)=1$, on $0\le t\le 8$. Its exact solution is the exponential decay $y(t)=e^{-t}$. Use step-doubling error estimation with the explicit Euler method: from the same numerical starting value $y_n$, form one coarse step $C_n$ and two half steps giving $F_n$:Accept $y_{n+1}=F_n$ when the absolute and relative error tolerances test $E_n=|{\hat{e}}_n|/s_n\le 1$ passes, with $s_n=10^{-5}+10^{-3}\max (|y_n|,|F_n|)$. There is no Richardson extrapolation in the accepted value. Both approximations have order of a numerical method one, while this local error estimator has step exponent $r=2$. Consequently the adaptive step-size controller uses a square root, not a cube root:Start with $h_0=1$. A rejected trial retains both $t_n$ and $y_n$, caps $q_n$ at one, and retries. A zero estimate uses $q_n=5$. Before every trial, cap $h_n$ at $8-t_n$. Stop with a failure if floating-point arithmetic gives $t_n+h_n=t_n$, or after 100000 trials. These rules specify the example completely; they are not a prescription for every ordinary differential equation solver.

The first two trials, $h=1$ and $h=0.2$, have $E\approx 247.525$ and $E\approx 9.90099$ and are rejected. The next trial has $h\approx 0.0572049$ and $E\approx 0.81$ and is accepted. For positive decreasing numerical values, $s_n=10^{-5}+10^{-3}{y}_n$, so the local error estimator givesWhen the relative term dominates, this suggests nearly constant accepted steps, approximately $2(0.9)\sqrt{10^{-3}}$. When $y_n$ falls below the crossover size $10^{-5}/10^{-3}=0.01$, the absolute term increasingly permits longer steps. The sharp final decrease in $h$ only lands on $t=8$; it does not indicate a sudden loss of smoothness.

The actual signed local defect from the same numerical input isHere $O$ is Big O notation, obtained from Taylor's theorem for the exponential function. Thus ${\hat{e}}_n$ agrees with the leading term as $h_n\to 0$, but is not identical to ${\ell }_n$. At $h=1$, $y_n=1$, the estimate is $0.25$ while the actual local defect is $e^{-1}-1/4\approx 0.117879$. Their closeness improves for the smaller accepted steps. Had we added ${\hat{e}}_n$ using Richardson extrapolation, the accepted formula would instead be $y_n(1-h_n+h_n^2/2)$, of order of a numerical method two; the same raw difference would still scale as $h_n^2$ and would not estimate that improved value's leading $h_n^3$ defect.

The global error $g_n=y_n-e^{-t_n}$ obeys an exact recurrence relation involving the actual local defects:Each previous global error contracts, yet the local defects keep adding. In the plotted run, 101 steps are accepted and two trials are rejected. Every accepted $E_n$ is at most approximately $0.81$, while the largest $|g_{n+1}|/s_n$ is approximately $38.6927$. At $t=8$, the numerical value is approximately $0.00027184218$, compared with the exact $e^{-8}\approx 0.00033546263$, giving absolute error approximately $6.36204\times 10^{-5}$. The right-hand plot deliberately scales actual global error by the same $s_n$ used for local acceptance; it is not a new global error tolerance test.

Tightening both absolute and relative error tolerances by a factor ten gives 315 accepted steps and endpoint absolute error approximately $2.14922\times 10^{-5}$; tightening both by a factor one hundred gives 991 accepted steps and endpoint absolute error approximately $6.94171\times 10^{-6}$. This is a convergence comparison for this particular initial value problem. The endpoint error decreases by about a factor three per tenfold tightening, rather than by ten: a first-order Euler method accumulates defects across many steps, while its local error estimator varies quadratically with $h$. A variable-step global ODE error bound needs propagation information as well as accepted local tests.

The absolute stability boundary remains a separate constraint. On the Dahlquist test equation $y^{\prime }=-\lambda y$ with $\lambda >0$, the accepted two-half-step formula has stability function $R(-\lambda h)=(1-\lambda h/2{)}^2$. It strictly damps perturbations only for $0<\lambda h<4$. Its local error estimator is ${\lambda }^2{h}^2{y}_n/4$, which can pass an absolute-tolerance test at tiny $y_n$ even outside this damping interval. In a stiff two-mode linear system, a fast mode can therefore constrain the step long after its visible amplitude has become small. An adaptive step-size controller does not turn an explicit Euler method into an A-stable method. Official solver documentation likewise distinguishes explicit Runge-Kutta methods for nonstiff problems from implicit Runge-Kutta methods and backward differentiation formulas for stiff equations.^[ref]

If a scaled local error estimator behaves as $E\propto h^r$ with approximately unchanged error scale, the displayed update targets a fixed local level. A safety factor $0<\eta <1$ aims below the acceptance boundary. For an embedded Runge-Kutta pair with order of a numerical method $p$ for the accepted output and $p-1$ for the companion formula, the difference generally has $r=p$; for the actual unscaled local truncation error of an order-$p$ formula, $r=p+1$. The exponent belongs to the quantity being controlled. Dividing the defect by $h$ changes that exponent. This follows by solving $E(h_{\mathrm{new}})\approx E(h)(h_{\mathrm{new}}/h{)}^r$ for the next step.^[ref]

Rejecting a trial means retrying from the same accepted time and state, with a smaller step. Advancing time or retaining the rejected output would change the initial value problem being approximated. Minimum and maximum change factors prevent extreme updates; a zero estimate requires finite capped growth. After acceptance, cap the next step at the remaining interval length. Detect a step that cannot advance time in floating-point arithmetic, and report failure rather than an endless retry. An adaptive step-size controller can use earlier estimates to smooth step sequences, but the memoryless rule here is sufficient to explain acceptance and rejection.^[ref]

For an initial value problem, convergence means that the largest global error on each fixed finite time interval tends to zero as the step size tends to zero and the starting errors vanish. For linear multistep methods, the Dahlquist equivalence theorem separates this requirement into consistency of a numerical method and zero-stability.

This initial value problem uses a square root and absolute value to define a continuous function on all of ${\mathbb{R}}^2$. For every finite waiting time $c\ge 0$, a forward solution of a differential equation on $t\ge 0$ isThe constant solution of a differential equation $y_{\infty }(t)=0$ is another choice. For example, at $t=3$, waits $c=0,1,2,3$ give values $9,4,1,0$. These are dimensionless mathematical trajectories; no measured physical waiting time is being supplied.

Verification uses the original ordinary differential equation, including the join. Before departure, both $y_c^{\prime }$ and $2\sqrt{|y_c|}$ vanish. For $t>c$, their common value is $2(t-c)$. For $c>0$, the difference quotient at the join isThus the derivative exists at $c$, equals the rate law there, and is continuous across the join. For $c=0$ the same computation gives the right derivative at the starting endpoint. For $c>0$, the left and right second derivatives at departure are $0$ and $2$: these forward solutions of a differential equation are continuously differentiable functions, without being twice differentiable at departure. The local two-sided statement of the Peano existence theorem is separate from this explicitly forward family.

The rate law tends to zero as $y\to 0$, but its secant slopes from zero are unbounded:No finite local Lipschitz constant exists at zero. At $\epsilon =10^{-2},10^{-4},10^{-6}$ the quotient is $20,200,2000$, even though the rates are shrinking to zero. Continuity controls the rate itself; Lipschitz continuity would control its change relative to the state displacement.

There is also a numerical lesson. The forward Euler method starting exactly at zero gives $Y_{k+1}=Y_k+2\Delta t\sqrt{|Y_k|}=0$ at every step. It reproduces the constant solution of a differential equation for every positive step size, while the moving trajectories remain equally valid. A computation producing one trajectory cannot establish uniqueness. Likewise, the separable differential equation calculation that divides by $\sqrt{|y|}$ only describes the moving pieces: checking the original equation recovers the zero pieces and their joins.

Evolution is deterministic when a specified initial condition and evolution rule select a unique future orbit. For a discrete dynamical system, applying the same map to the same state gives the same next state. For an ordinary differential equation, uniqueness of the initial value problem is required; merely writing a differential equation does not guarantee it.

Determinism permits both asymptotic stability and sensitive dependence on initial conditions. It does not by itself imply reliable prediction from approximate measurements, and sensitivity does not introduce a random choice into the evolution rule.

A direction field for $y^{\prime }=f(t,y)$ places a short segment with slope $f(t,y)$ at each sampled point where $f$ is defined. A function graph of a solution of a differential equation follows these segments: its tangent line has exactly the prescribed slope. Segment length is a drawing choice, rather than a speed; arrows, if used, point toward increasing $t$. An isocline collects points with one common slope.

For Newton's law of cooling, $y^{\prime }=-0.1(y-20)$ with time in minutes and temperature in degrees Celsius, the slopes at $y=80,40,20,10$ are respectively $-6,-2,0,1$ degrees Celsius per minute. Each horizontal row is an isocline, because this autonomous differential equation has no explicit time dependence. Above $20$ the solutions of a differential equation descend, below $20$ they rise, and $y=20$ is an equilibrium point.

The highlighted initial condition $y(0)=80$ selects one curve from the family. The picture reveals approach to the equilibrium point before any algebra is done. Under the hypotheses of the Picard-Lindelöf theorem, distinct function graphs of solutions cannot cross or even meet in the time-state plane. A direction field alone does not guarantee uniqueness: several curves can share a tangent line at the same point when those hypotheses fail. See the delayed-start example under initial value problem.

This geometric interpretation is developed in Direction Fields, Isoclines, Integral Curves.

Take a reservoir with constant horizontal cross-sectional area $A$, one small orifice of area $a$, no inflow, and constant stipulated discharge coefficient $C_d>0$. Let the initial head be $H>0$. For constant mass density $\rho$, the stored mass above the outlet level is $\rho Ah$ and outgoing mass flow rate is $\rho Q$. Mass conservation cancels $\rho$ and gives the initial value problemThe downward speed of the free surface is $Q/A$, smaller than the ideal jet speed by the factor $C_{d}a/A$. Thus a small opening relative to the reservoir supports the negligible-surface-velocity assumption.

For $h>0$, separation of variables or the chain rule gives $d\sqrt{h}/dt=-k/2$. Integrating yieldsThe model reaches zero head in finite time even though volumetric flow rate tends to zero. Extending the squared expression beyond $T$ would make head grow again and violate the original ordinary differential equation. The zero-head continuation enforces no further outflow. For an elevated side opening, zero head means the surface has reached that opening, leaving fluid below it; $T$ is not the time to remove every drop from such a tank.

For an original synthetic example, stipulateThe initial volume above outlet level is $AH=0.0250{\mathrm{m}}^3=25.0\mathrm{L}$. The initial ideal jet speed is $2.215\mathrm{m}{\mathrm{s}}^{-1}$, the initial volumetric flow rate is $Q_0=0.1373\mathrm{L}{\mathrm{s}}^{-1}$, and the predicted zero-head time is $T=364.1\mathrm{s}$. These are consequences of stipulated parameters, not experimental measurements. At $t=T/2=182.1\mathrm{s}$, head is $H/4=0.0625\mathrm{m}$, stored volume is $6.25\mathrm{L}$, and volumetric flow rate is $Q_0/2=0.06866\mathrm{L}{\mathrm{s}}^{-1}$. Half the initial volume has already left at $t=T(1-1/\sqrt{2})=106.7\mathrm{s}$.

An independent mass conservation check integrates the outgoing volumetric flow rate:Multiply both sides by mass density to obtain discharged mass. A dimensional analysis check gives $[k]={\mathrm{m}}^{1/2}{\mathrm{s}}^{-1}$, $[Q]={\mathrm{m}}^3{\mathrm{s}}^{-1}$ and $[T]=\mathrm{s}$. The triangular area under $Q(t)$ explains why the initial-rate estimate $AH/Q_0$ is only half the actual model time. Doubling $A$ doubles $T$; doubling $a$ or $C_d$ halves $T$, within the same assumptions.

This finite-time result is an idealized prediction. The constant $C_d$ and quasi-steady approximation need not survive very small head or the initial startup. Appreciable viscosity in a narrow outlet tube requires a different pressure balance, as in Poiseuille flow; appreciable surface tension can prevent zero head altogether.

For $x>0$ on a domain where the tangent is finite, put $v=y/x$. The homogeneous first-order differential equation becomes $x{v}^{\prime }=-\tan v$. Direct differentiation gives $(x\sin v{)}^{\prime }=0$, hence $x\sin (y/x)=C$, without dividing by $\sin v$ and losing constant-ratio solutions. If $C=x_0>0$ and $y/x\to \pi /2$ as $x\downarrow x_0$, two branches are $y=x\arcsin (x_0/x)$ and $y=x[\pi -\arcsin (x_0/x)]$, defined for $x>x_0$. They approach the same finite endpoint but have unbounded derivatives. The original tangent is undefined at the endpoint, so neither is a classical solution through that point. An implicit relation extending continuously to an initial point does not by itself solve an initial value problem there.

For an incoming timelike Schwarzschild geodesic with asymptotic Lorentz factor $k\ge 1$, the critical unstable circular orbit obeysAt $k=1$, $r_c=4\mu$ and $|h_c|=4\mu c$. For $k>1$, the critical radius lies between $3\mu$ and $4\mu$ and $|h_c|>4\mu c$. The marginally bound angular-momentum bound remains sufficient at larger energies but is no longer the exact threshold. Origin at infinity and an inward initial branch are essential; an orbit already inside a barrier is a different initial value problem.

Consider identical harmonic oscillators with unit mass, natural angular frequency ${\omega }_0>0$, and positive coupling $\kappa \ll {\omega }_0^2$:An initial impulse on the first gives $x_1(0)=x_2(0)={\dot{x}}_2(0)=0$ and ${\dot{x}}_1(0)=v$. The normal mode coordinates $u=x_1+x_2$ and $w=x_1-x_2$ have angular frequencies ${\omega }_0$ and ${\omega }_a=\sqrt{{\omega }_0^2+2\kappa }$. Solving their initial value problems provesLet $\Delta ={\omega }_a-{\omega }_0\sim \kappa /{\omega }_0$. To leading order in $\kappa /{\omega }_0^2$, sum-to-product identities give the displacement amplitude envelopes $(|v|/{\omega }_0)|\cos (\Delta t/2)|$ and $(|v|/{\omega }_0)|\sin (\Delta t/2)|$. Thus the first nearly complete energy transfer occurs at $\pi /\Delta$, and localization returns at $2\pi /\Delta$. More precisely, at $t=\pi /\Delta$, the first oscillator has exactly zero velocity and residual displacement $\frac{v}{2}(1/{\omega }_0-1/{\omega }_a)\sin ({\omega }_{0}t)$, of order $|v|\kappa /{\omega }_0^3$.

Here transferred energy means kinetic energy plus the uncoupled quadratic potential energy; the coupling potential energy is relatively small. The fast phase need not make the receiving oscillator pass through zero displacement at the transfer time. Expanding $\Delta t$ as small is invalid on the transfer scale $t=O(1/\kappa )$. This lemma complements absolute and signed beat envelopes by keeping that slow phase.

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.

The Euler method approximates an initial value problem $y^{\prime }=f(t,y)$, $y(t_0)=y_0$, by evaluating the derivative at the current point and moving along that tangent line for a time $h$. With $t_n=t_0+nh$, its displayed update needs one evaluation of $f$ per step. This simple rule turns local slope information into a whole approximate trajectory. MIT's course notes derive it by dropping the second-order Taylor remainder.^[ref]

Take $y^{\prime }=-2y$, $y(0)=1$, whose exact solution is the exponential function $y(t)=e^{-2t}$. The Euler method givesAt $h=0.2$, the first values are $1,0.6,0.36,0.216$. At $t=1$, five steps give $0.6^5=0.07776$, compared with $e^{-2}\approx 0.135335$. The result decays but is too small: absolute stability alone does not ensure accuracy. At $h=1.2$, the first values are $1,-1.4,1.96,-2.744$, and at $t=6$ five steps give $-5.37824$, while the exact solution is about $0.00000614421$.

For the Dahlquist test equation $y^{\prime }=\lambda y$, the Euler method has stability function $R(z)=1+z$ with $z=h\lambda$. Its absolute stability region is $|1+z|\le 1$, the closed disk centred at $-1$ with radius one. For the present real decay problem and $h>0$, bounded numerical values require $0<h\le 1$; strict decay requires $0<h<1$. At $h=1$ the values alternate $1,-1,1,\dots$ without decay. These conditions follow directly from the amplification factor $1-2h$, matching MIT's criterion for $y^{\prime }=-ay$.^[ref]

The Euler method is first order: with bounded second derivative and suitable Lipschitz continuity, its global error is $O(h)$ on a fixed interval. In this example the logarithm gives $\log (1-2h)=-2h-2h^2+O(h^3)$ for small positive $h$. At $t=nh$, this means $y_n=e^{-2t}(1-2th+O(h^2))$. It explains both the downward bias and why halving the step roughly halves the error once $h$ is sufficiently small.

Under the Picard-Lindelöf theorem's hypotheses on $U\subset \mathbb{R}\times {\mathbb{R}}^d$, extend a solution of a differential equation as far as possible, to a maximal interval $(T_{-},T_{+})$. If $T_{+}<\infty$, the time-state trajectory $(t,y(t))$ must eventually leave every compact subset of $U$. Equivalently, confinement near $T_{+}$ to one compact subset of $U$ permits continuation past $T_{+}$. See Teschl, Section 2.6, Lemma 2.14 and Corollaries 2.15–2.16, pp. 52–53.

The domain restriction matters even for a bounded trajectory. Restrict $y^{\prime }=1$, $y(0)=0$, to $U=\mathbb{R}\times (-1,1)$. Then $y(t)=t$ exists in this domain only for $-1<t<1$; at $t=1/2,3/4,9/10$ its values are $1/2,3/4,9/10$. As $t\uparrow 1$, it approaches the excluded boundary point $(1,1)$ while remaining bounded. Its trajectory cannot fit into any compact subset of $U$ near that endpoint. This is escape from a finite domain, without finite-time blowup in norm. Enlarging the domain to $\mathbb{R}\times \mathbb{R}$ gives the same formula $y(t)=t$ for every real $t$.

Conversely, $y^{\prime }=y$, $y(0)=1$, stays in the compact set $[0,T]\times [1,e^T]$ for every fixed $T>0$. Its growing exponential function is unbounded over the entire future, but has no finite-time blowup. At $T=2$, its largest value on $[0,2]$ is only $e^2\approx 7.38905610$. The relevant question is boundedness before each finite endpoint, rather than one bound valid forever. The tangent function in the diagram under the Picard-Lindelöf theorem illustrates the other possibility: its norm escapes every bound before $T=\pi /2$.

This criterion is useful because an estimate can secure continuation without an explicit formula. In a finite-dimensional vector space, a bounded closed ball is compact; a bounded state trajectory therefore supplies the needed compact set when its time-state closure is inside the domain. That final proviso prevents mistaking a domain boundary for finite-time blowup.

A common sufficient estimate comes from global Lipschitz continuity in the state. Assume $U=I\times {\mathbb{R}}^d$ for an open time interval $I$, with $f$ continuous and a Lipschitz constant $L\ge 0$ valid for all states on each compact time interval under consideration. For $[t_0,T]\subset I$, let $B={\sup }_{t_0\le s\le T}\Vert f(s,0)\Vert$, finite by continuity. Then $\Vert f(t,y)\Vert \le B+L\Vert y\Vert$, and the Gronwall inequality gives, for $t_0\le t\le T$ wherever the solution already exists,For instance, $y^{\prime }=y+1$, $y(0)=0$, has $L=B=1$ and the exact solution $y(t)=e^t-1$, attaining this estimate; at $t=1$ the value is $e-1\approx 1.71828183$. Finite bounds of this kind give continuation throughout $I$, with the analogous estimate backwards in time. When $I=\mathbb{R}$, this gives existence for all real times. If $I$ has a finite endpoint, it gives no continuation beyond that prescribed domain. The underlying linear-growth criterion is ^[ref].

A brief explanation of the continuation mechanism suffices: within the protecting compact set, the continuous function $f$ is bounded, say by $M$. The integral equation yields $\Vert y(t)-y(s)\Vert \le M|t-s|$ near the finite endpoint, so $y(t)$ has a limit $y_{\ast }$. The protecting compact set puts $(T_{+},y_{\ast })$ inside $U$. Restarting the initial value problem there and using local uniqueness joins the trajectories.

For a system of linear differential equations $u^{\prime }=A(t)u$, a fundamental matrix has columns forming a basis of its solution space. It solves $X^{\prime }=A(t)X$ and is invertible at every time in its interval. Every solution is $u(t)=X(t)c$ for a constant vector $c$: differentiating $X^{-1}u$ gives zero. Invertibility at one time implies invertibility throughout the interval by uniqueness of the initial value problem. This notion is distinct from a Markov-chain fundamental matrix or the projective-geometric matrix used in stereo vision.

Hyperbolic functions combine the exponential function to describe a hyperbola much as trigonometric functions describe a circle. Their two basic members are the hyperbolic cosine and hyperbolic sine, defined for every real number $t$ byTheir derivatives exchange the two functions:These formulas follow by applying differentiation to the exponential function definitions. ^[ref]

Adding and subtracting gives $\cosh t+\sinh t=e^t$ and $\cosh t-\sinh t=e^{-t}$. Multiplying these equalities proves the hyperbolic Pythagorean identity:Consequently the parametric equation $(x,y)=(\cosh t,\sinh t)$ traverses the entire right branch of the unit hyperbola $x^2-y^2=1$, where $x\ge 1$. For example, $t=0$ gives $(1,0)$, while $t=\pm 1$ gives approximately $(1.543,\pm 1.175)$. Compare $(\cos \theta ,\sin \theta )$ on the unit circle, where the squared coordinates add to one. The minus sign changes the geometry from a closed circle to an unbounded hyperbola.

The other members are ratios and reciprocals: hyperbolic tangent is $\tanh t=\sinh t/\cosh t$, hyperbolic cotangent is $\coth t=\cosh t/\sinh t$, and hyperbolic secant is $\mathrm{sech}t=1/\cosh t$; denominators must be nonzero. The hyperbolic sine and hyperbolic cosine both solve the second-order linear differential equation $f^{\prime \prime }=f$. In fact $f(t)=A\cosh t+B\sinh t$ has $f(0)=A$ and $f^{\prime }(0)=B$, making this basis convenient for an initial value problem. Their geometry also makes hyperbolic substitution useful for integrals containing $\sqrt{1+u^2}$, and the hyperbolic cosine gives the physical catenary of a hanging chain.

An implicit ordinary differential equation relates the derivative to the independent and dependent variables without requiring it to be explicitly isolated. For example $y{y}^{\prime }=1$ is meaningful at every pair $(x,y)$ as an equation, but no differentiable solution of a differential equation can pass through $y=0$: its left side would vanish there. Dividing by $y$ gives $y^{\prime }=1/y$ only on the domain $y\ne 0$. Conversely multiplying an explicit equation by a vanishing denominator can introduce points where the original equation was undefined. A solution of the multiplied equation at such a point need not solve the original initial value problem.

An implicit solution describes a solution of a differential equation through a relation $F(x,y)=C$, instead of isolating $y$ in an explicit formula. The relation describes a level set; an actual scalar solution of a differential equation must be a differentiable function $y(x)$ on an interval whose function graph lies in the equation's domain. A whole level set need not be one such function graph.

For example, considerImplicit differentiation of $x^2+y^2=C$ gives $2x+2y{y}^{\prime }=0$. For $C>0$, the circle contains two solutions of a differential equation,Both satisfy the equation on the stated open interval, and their function graphs stay in $D$. An initial condition $(x_0,y_0)\in D$ fixes $C=x_0^2+y_0^2$ and selects the sign of $y_0$. In particular, $y(0)=2$ gives $C=4$ and the upper branch on $(-2,2)$. The constant alone does not select the branch.^[ref]

The precise local criterion uses the implicit function theorem. Let $F$ be continuously differentiable on an open subset of the plane, with $F(x_0,y_0)=C$. If its partial derivative $F_y(x_0,y_0)\ne 0$, the nearby level set is a unique continuously differentiable function graph $y=g(x)$ through that point. The chain rule givesThis proves that the local function graph solves $y^{\prime }=-F_x/F_y$ there. For a different prescribed equation $y^{\prime }=f(x,y)$, one must also check $F_x+F_{y}f=0$ along that branch and that it stays in the domain of $f$. A level set by itself does not verify an arbitrary ordinary differential equation.^[ref]

For $F=x^2+y^2$, the upper branch through $(1,\sqrt{3})$ has $F_y=2\sqrt{3}$ and tangent line slope $-1/\sqrt{3}$. At $(2,0)$, however, $F_y=0$ and $F_x=4$. The circle remains a smooth regular level set with a vertical tangent line: it can be written locally as $x=\sqrt{4-y^2}$, but not as a single function graph $y=g(x)$ over an open interval containing $2$. Moreover, $-x/y$ is undefined there, so the scalar solution of a differential equation cannot include this point in its domain. Changing the parameter draws the curve through the point but does not extend that scalar initial value problem.

When both partial derivatives vanish, even the geometry needs separate inspection. The level $x^2+y^2=0$ is an isolated point, not a function graph on any nonempty open interval. The zero level set of $G(x,y)=y^2-x^2$ instead consists of the two lines $y=x$ and $y=-x$ crossing at the origin, where $\nabla G=0$. Each line satisfies the implicit ordinary differential equation $y{y}^{\prime }=x$, including at the origin; the explicit form $y^{\prime }=x/y$ excludes the origin. Thus algebraically dividing by a vanishing coefficient changes the domain and can change which solutions of a differential equation are admissible. Vanishing $F_y$ alone is only failure of a sufficient test: $F(x,y)=y^3$ has the smooth zero function graph $y=0$ despite $F_y=0$ everywhere on it.

Implicit solutions are especially useful for exact first-order ordinary differential equations and for finding a first integral of an ordinary differential equation: a conserved relation can reveal allowed states without an explicit formula. For the normalized harmonic oscillator written as $\dot{x}=y$, $\dot{y}=-x$, the energy $H=(x^2+y^2)/2$ has $\dot{H}=x\dot{x}+y\dot{y}=0$. Its positive level sets are the same circles. Where $y\ne 0$, eliminating time gives $dy/dx=-x/y$ and the local implicit solution above. The harmonic oscillator's parametric curve continues through $y=0$ because its vector field governing evolution in time is defined there; treating $x$ as the independent variable fails at that turn. The origin is an equilibrium point of the harmonic oscillator, although it is excluded from the scalar equation $y^{\prime }=-x/y$.

An initial value problem combines an ordinary differential equation with an initial condition, for example $y^{\prime }=f(t,y)$ and $y(t_0)=y_0$. The task is to find a solution of a differential equation on an interval containing $t_0$ which has that starting value. For a continuous linear differential equation $u^{\prime }=A(t)u$, every prescribed finite vector $u_0$ determines a unique solution on the coefficient interval. A boundary value problem can instead prescribe values at distinct endpoints.

For the cooling initial value problem $y^{\prime }=-k(y-T_a)$, $y(t_0)=y_0$, where $k>0$ and $T_a$ are constant, the unique solution isTake $T_a=20$ degrees Celsius, $y_0=80$ degrees Celsius, $t_0=0$, and $k=0.1$ per minute. Then $y(10)=20+60e^{-1}\approx 42.07$ degrees Celsius. Reaching $30$ degrees Celsius requires $60e^{-0.1t}=10$, so $t=10\log 6\approx 17.92$ minutes. The time constant $1/k=10$ minutes reduces the excess temperature by a factor $e^{-1}$; it does not reduce the absolute temperature by that factor.

Uniqueness needs assumptions. If $f$ is a continuous function on an open neighborhood of $(t_0,y_0)$ and satisfies a local Lipschitz condition in $y$, uniformly in $t$, the Picard-Lindelöf theorem gives a unique local solution. Concretely, a finite $L$ bounds $|f(t,u)-f(t,v)|\le L|u-v|$ in a sufficiently small rectangle. Cooling satisfies this with $L=k$. Local existence does not promise existence for all future time: $y^{\prime }=y^2$, $y(0)=1$, has $y(t)=1/(1-t)$ on $(-\infty ,1)$ and finite-time blowup at $t=1$.

Continuity alone is insufficient for uniqueness. Consider the dimensionless initial value problem involving a square root and an absolute value:Besides $y(t)\equiv 0$, each delay $a\ge 0$ gives a solution for $t\ge 0$:At $t=4$, delays $a=0,1,2$ give $y=4,9/4,1$, despite identical starting data. These are valid solutions of a differential equation: for $a>0$, both one-sided derivatives at the join are zero; for $a=0$, the right derivative at the starting time is zero. After the join, $y_a^{\prime }=(t-a)/2=\sqrt{|y_a|}$ because $t>a$. Squaring alone would also suggest a parabola before the join, but its negative derivative there would violate the equation. Since $y^{\prime }\ge 0$, a forward solution starting from zero stays nonnegative.

The right-hand side is a continuous function but fails the Lipschitz condition at zero: for $h>0$, $|\sqrt{|h|}-0|/|h|=1/\sqrt{h}$ is unbounded. Thus the Peano existence theorem still supplies local existence, while the Picard-Lindelöf theorem does not supply uniqueness. Related waiting-time examples and a sufficient uniqueness condition appear in MIT 18.034, Lecture 2: Fundamental Principles.

A brief proof outline explains the mechanism: rewrite the initial value problem as the integral equation $y(t)=y_0+{\int }_{t_0}^{t}f(s,y(s))ds$. On a short enough interval, the Lipschitz condition makes the right-hand side a contraction mapping on a closed set of continuous functions near $y_0$. The contraction mapping theorem gives a unique fixed point; the fundamental theorem of calculus turns it back into a solution.

An integrating factor is a nonzero multiplier that makes a differential equation directly integrable. For a first-order linear differential equation, it gathers the unknown function and its derivative into one product rule expression. Suppose $a$ and $b$ are continuous functions on an interval containing $t_0$, and consider the initial value problemDefine $A(t)={\int }_{t_0}^{t}a(s)ds$ and $\mu (t)=e^{A(t)}$. The exact solution of a differential equation throughout that interval isThe definite integral fixes $\mu (t_0)=1$, so the initial condition appears without an extra constant to determine. If the equation originally has a coefficient multiplying $y^{\prime }$, first divide by that coefficient on an interval where it is nonzero; a zero of that coefficient needs separate attention. See OpenStax's first-order linear equations chapter for this normalization and the integrating factor procedure.

For a concrete balance, solve $y^{\prime }+2y=3$ with $y(0)=0$. Here $\mu =e^{2t}$, so multiplying the entire ordinary differential equation gives $(e^{2t}y{)}^{\prime }=3e^{2t}$. Accumulating from zero to $t$ givesThe computed values are $y(0.5)\approx 0.94818$, $y(1)\approx 1.29700$ and $y(2)\approx 1.47253$. Its derivative is $3e^{-2t}$, which verifies $y^{\prime }+2y=3$ exactly. Interpreting $3$ as a steady input and $2y$ as proportional removal explains the shape: initially all input increases $y$, then removal catches up. The constant particular solution $y_p=3/2$ balances input and removal, and the homogeneous solution $y_h=-\frac{3}{2}{e}^{-2t}$ is the transient response that corrects it to the initial condition. Because this autonomous differential equation has a time-independent right-hand side, $3/2$ is also an equilibrium point.

The first panel adds the particular solution and homogeneous solution; the second shows the same balance in rates, $3=y^{\prime }+2y$. This links the algebra to a physical accounting rule. Newton's law of cooling has this form with $a=k$ and $b=k{T}_a$, while an electrical circuit can have a varying input. In both cases the integrating factor systematically handles the input rather than requiring a guess for a particular solution.

Variable coefficients work just as directly. For $y^{\prime }+2ty=t$, $y(0)=0$, the integrating factor is $e^{t^2}$, andThus $y(0.5)\approx 0.11060$, $y(1)\approx 0.31606$ and $y(2)\approx 0.49084$. The derivative $y^{\prime }=t{e}^{-t^2}$ verifies $y^{\prime }+2ty=t$, including at zero, where both coefficient and input vanish. No division by $t$ is necessary. The general solution is $y=1/2+C{e}^{-t^2}$: $1/2$ is a constant particular solution, and $C{e}^{-t^2}$ solves the associated homogeneous linear differential equation. Although this special forcing term admits a constant trajectory, this is a nonautonomous differential equation: it depends explicitly on time. Merely setting its right-hand side to zero at one instant would not establish a constant solution of a differential equation; at $t=0$ every state has zero slope.

The direction field shows the selected trajectory leaving its initially horizontal tangent line. For $t\ge 0$, the difference between any two solutions of a differential equation is multiplied by $e^{-t^2}$, rather than $e^{-2t}$ in the constant-rate example. This dependence on elapsed history is clearer if the exact formula is rewritten asFor $t\ge t_0$, each contribution from the input at $s$ is weighted by the accumulated coefficient between $s$ and $t$. The homogeneous solution remembers the initial condition; the particular solution given by the second term accumulates the forcing term with zero initial value. This decomposition differs from choosing a constant particular solution, but gives the same total. Initial-data memory tends to zero when ${\int }_{t_0}^{t}a\to +\infty$; a negative accumulated coefficient instead amplifies it. The reciprocal relation between integrating factor and homogeneous solution is developed in MIT's notes on solutions to linear first-order equations.

Why does the trick work? The fundamental theorem of calculus and chain rule give ${\mu }^{\prime }=a\mu$. Consequently the product rule givesIntegrating this identity and dividing by the positive exponential function $\mu$ produces the stated formula. Any nonzero constant multiple of $\mu$ works equally well; normalization is a convenience. This is the beauty of the method: a sum involving an unknown function becomes one exact derivative, leaving only an integral of known functions. See Oregon State's product-rule derivation. The related integrating factor for a differential one-form extends the idea of making an expression exact beyond this linear setting.