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]

The Null Raychaudhuri equation is the light-ray version of the Raychaudhuri equation associated with Amal Kumar Raychaudhuri. His original work concerned relativistic cosmology.^[ref] For the normal null geodesic congruence before its first focal point, null twist vanishes. With an affine parameter $\lambda$, the four-dimensional equation reduces to the displayed ordinary differential equation. The null shear contraction is nonnegative on the transverse screen space, and the null convergence condition makes the other part of $Q$ nonnegative. Consequently ${\theta }^{\prime }\le -{\theta }^2/2$.

For $\theta (0)={\theta }_0<0$, put $w=-1/\theta >0$. ThenA regular negative null expansion therefore cannot persist beyond $\lambda =2/|{\theta }_0|$ if the null geodesic extends that far. The null focusing theorem turns this estimate into a focal-point bound. The number refers to an affine parameter, not the proper time of a light ray, which is zero.

For an original soluble model, prescribe constant $Q$ along one twist-free beam; this does not specify an entire solution of the Einstein field equations. Introduce a positive affine scale $L$ and the dimensionless variablesFor an infinitesimal area $A$ define $a=\sqrt{A/A_0}$ before the first zero. Since $y=2a^{\prime }/a$, where primes now mean $d/du$, the same model becomesTwo independently checkable solutions areThe nonnegative extra focusing term moves the first zero earlier. Direct substitution into both ordinary differential equations verifies the solutions; numerical integration of the nonlinear equation provides a separate comparison with the linear beam equation.

A focal point makes this beam description degenerate, but individual null geodesics may continue smoothly through it. Even converging light in Minkowski spacetime can form a caustic. Converting local focusing into null geodesic incompleteness requires the global step below.

This second-order ordinary differential equation has the independent entire Airy functions $\mathrm{Ai}$ and $\mathrm{Bi}$ as a standard basis. It is distinct from the dispersive PDE in Airy equation. Its general solution is $a_1\mathrm{Ai}(Z)+a_2\mathrm{Bi}(Z)$.

Take the following bounded analytic example for a dimensionless positive coupling $g$, with fixed $ϵ>0$ and $g(0)=g_0>0$:Increasing $\ell$ means increasing observation length $L$, hence motion toward the infrared fixed point. This convention matters: if the inverse-length scale is $\mu ={\mu }_0{e}^{-\ell }$, the same trajectory has renormalization-group beta function ${\beta }_{\mu }(g)=\mu dg/d\mu =-ϵg+g^2$. The two derivatives have opposite signs.

To solve the equation, set $v=1/g$. It obeys $dv/d\ell =-ϵv+1$, a linear ordinary differential equation. ThereforeFor $\ell \ge 0$ its denominator is positive, and $g$ stays between $g_0$ and $ϵ$. Thus the trajectory is bounded on this domain. The renormalization-group fixed points are $g=0$ and $g_{\ast }=ϵ$. The derivative of the length-flow right-hand side is $ϵ-2g$: zero is repulsive from positive $g$, whereas $g_{\ast }$ has eigenvalue $-ϵ$ and attracts positive trajectories. Linearization gives $g-g_{\ast }\propto e^{-ϵ\ell }=b^{-ϵ}$ near $g_{\ast }$. The inverse-length beta function has positive derivative at the same infrared-attractive point.

For the chosen $ϵ=1/5$, take two starting models with $g_0=1/20$ and $g_0=3/5$. At $b=1024$, $e^{-ϵ\ell }=1/4$, so their couplings are respectively $4/35$ and $6/25$. Both move toward $1/5$, from opposite sides. They are still visibly different: slow decay of an irrelevant operator can leave substantial corrections at accessible lengths. The case $g_0=0$ stays exactly zero and is outside this positive-coupling basin.

An axisymmetric Stokes streamfunction disturbance $\varphi (r)e^{ik(z-ct)}$ of an axial basic flow obeys the displayed ordinary differential equation. Regularity requires $\varphi =O(r^2)$ at the axis, and an impermeable wall requires $\varphi (a)=0$ for nonzero $k$. Multiplication by ${\varphi }^{\ast }/r$ and integration by parts shows that a growing real-wavenumber mode requires $(U^{\prime }/r{)}^{\prime }$ to change sign; if it vanishes identically, the positive energy integral rules out such a mode. The condition concerns this cylindrical vorticity gradient rather than $U^{\prime \prime }$ alone.

The Beer-Lambert law relates absorbance $A$ to an absorbing species' concentration $c$ and path length $l$ by $A=\epsilon lc$ in a homogeneous medium using sufficiently narrow-band light. The coefficient $\epsilon$ depends on the species, wavelength and chemical solution conditions. This makes Ellman's reagent useful for quantitative thiol assays once its released colored species has been calibrated. ^[ref]

Here $c$ is molar concentration, $l$ is the geometric length through the sample, and $\epsilon$ is the molar absorption coefficient. The ideal beam has parallel light rays, a single wavelength, and a medium with uniform composition and isotropy; finite spectral width must be small relative to features in the optical spectrum.^[ref]

For a spectrophotometer, distinguish the source's output from the reference that reaches the photodetector. Let $P_0$ be the received radiant power through a reference blank and $P$ the received radiant power through the sample, after subtracting the photodetector's dark response. Both refer to the same selected wavelength band and collection geometry. Define the blank-normalized transmittance and absorbance byThe reference blank includes the solvent and other assay ingredients except the absorbing target. Matching cuvettes and conditions makes their common losses divide out. Thus this $T$ is relative to the reference blank; it is not the absolute transmittance of the entire cuvette and its contents. Both radiant powers can be expressed in watts, while their ratio, $T$, and its logarithm, $A$, are dimensionless. An instrument displaying $40\%$ transmission means $T=0.40$, not $T=40$.

Why take a logarithm? Consider two consecutive identical absorbing layers. Each transmits a fraction $t$ of the radiant power entering it, so together they transmit $t^2$. Their absorbances add because $-{\log }_{10}(t^2)=2(-{\log }_{10}t)$. Doubling the sample length doubles absorbance while squaring transmittance. This turns multiplicative loss into an additive scale and makes a straight line useful for inferring molar concentration.

A short model derivation makes the assumption explicit. Suppose a thin layer of thickness $dx$ removes the fraction $kcdx$ of its incoming radiant power, where $k$ is a fixed constant of proportionality and $c$ is a uniform molar concentration. The corresponding ordinary differential equation, initial condition, and solution areTaking the base-ten logarithm gives $A=kcl/\ln 10$, so $\epsilon =k/\ln 10$. The exponential decay follows from equal layers removing equal fractions of what remains. Proportionality to molar concentration additionally assumes that the absorbing chemical species retains the same response as its molar concentration changes; the geometric argument alone does not establish that.

A Bernoulli differential equation is a first-order ordinary differential equation of the displayed form, with $P$ and $Q$ continuous functions on an interval and a fixed real exponent $n$. A nonlinear power can hide a first-order linear differential equation: on a suitable branch, the Bernoulli power substitution $z=y^{1-n}$ removes it exactly. The transformed solution of a differential equation can be completely regular even when recovering $y$ produces finite-time blowup.

For example, the reciprocal Bernoulli equation with finite-time blowup $y^{\prime }+y=e^x{y}^2$ becomes $z^{\prime }-z=-e^x$ under $z=1/y$. The initial conditions $y(0)=1/2$ and $y(0)=-1$ select a positive trajectory blowing up at $x=2$ and a negative trajectory tending to zero as $x\to \infty$, respectively. The integrating factor exposes both behaviors with the same short calculation.

The cases $n=0$ and $n=1$ are already first-order linear differential equations: respectively $y^{\prime }+Py=Q$, with the zeroth-power term interpreted as the constant $1$, and $y^{\prime }+(P-Q)y=0$. There is no need for a nonlinear change of variable; for $n=1$, $y^{1-n}=1$ would discard the unknown entirely. For other exponents, the substitution must respect the domain of exponentiation and may lose the zero solution of a differential equation.^[ref]

The biomass concentration is the mass of biological material per unit culture volume; here it is measured as dry cell mass in $\mathrm{g}{\mathrm{L}}^{-1}$. A biomass concentration differs from a cell count per unit volume because cell sizes can change. In the ideal chemostat, biomass leaves the vessel at rate $Qx$, giving a removal term $-Dx$ in the ordinary differential equation for $x$.

A bounded function on $[0,\infty )$ has a Laplace transform satisfying the displayed bound, by integration of $|y(t)e^{-pt}|$. Consequently its Laplace transform cannot have a pole of order two or higher at $p=0$. This can eliminate the arbitrary constant introduced when multiplication by $t$ turns an ordinary differential equation into a first-order differential equation for its Laplace transform.

For example, transforming $t{y}^{\prime \prime }+(t-1)y^{\prime }-y=0$ with $y(0)=1$ gives $p(p+1)Y^{\prime }+(3p+2)Y=2$. An integrating factor $p^2(p+1)$ yieldsThe bound forces $C=0$, selecting $y=e^{-t}$. This selection can be justified without assuming unproved endpoint regularity for every solution: $e^{-t}$ and $t-1$ have nonzero Wronskian $t{e}^{-t}$ on $(0,\infty )$, so they form a fundamental system of solutions. Boundedness removes $t-1$, and the limit at zero fixes the remaining coefficient. The endpoint limit by itself does not select a unique solution.

A box model represents a system by finitely many internally uniform compartments. Balances of transported material and local sources give coupled ordinary differential equations. For a constituent's concentration $c_i$ in fixed box volume $V_i$, the stored amount is $V_i{c}_i$; changing volumes requires the balance of that product. A transport-only closed box model preserves total inventory, even when its transport matrix is nonsymmetric.

The butterfly effect is the intuitive idea that an arbitrarily small change in an initial condition can eventually change the observed future substantially. In dynamical systems, sensitive dependence makes this precise for a map $F:X\to X$ on a metric space $(X,d)$: there is one separation scale $\delta >0$ such thatThe same $\delta$ works at every point, whereas $y$ and the separation time $n$ may depend on the point and the requested closeness. For a flow map, replace the iteration of a map $F^n$ by ${\Phi }_t$ with some $t\ge 0$. The state space matters: restricting to a single fixed point removes the alternative nearby states needed for sensitivity.

Sensitivity is compatible with deterministic evolution. Given an exact initial condition, there is still one orbit; uncertainty concerns which exact state a measurement represents. It also allows accurate short-term prediction. For a continuous function $F$, each fixed finite iterated function $F^n$ is continuous, so sufficiently close initial states remain close at any specified finite step. What fails is a single initial accuracy that controls every future step.

The logistic map illustrates both sides. At $r=2.8$, nearby orbits in the basin of attraction of its attracting fixed point approach one another. At $r=4$, arbitrarily close states can eventually separate by a macroscopic fraction of $[0,1]$. The definition asserts the existence of separating neighbours, not that every nearby pair must separate: $x$ and $1-x$ have exactly the same next state under the logistic map and hence the same subsequent orbit. Nor must separation increase at every step.

A positive Lyapunov exponent describes average exponential growth of an infinitesimal displacement along an orbit, when the relevant limit exists. It is related to sensitivity but is not a substitute for its uniform definition. For a differentiable interval map, the chain rule givesAt $r=4$, the logistic map has $\lambda =\log 2$ for Lebesgue measure-almost every initial condition, while exceptional orbits can have different exponents. Finite-separation growth is only approximately exponential while the linear approximation remains useful. In a bounded state space, a separation cannot grow exponentially forever.

Sensitivity alone is also weaker than Devaney chaos, which additionally requires topological transitivity and dense periodic points. For example, $F(x)=2x$ on the real line has sensitive dependence and a positive Lyapunov exponent, but its only periodic point is zero, so it is not Devaney-chaotic. This distinguishes growth of uncertainty from the richer organization of recurrent orbits.

Edward N. Lorenz's "Deterministic Nonperiodic Flow" (1963) studied a thermal convection model with three ordinary differential equations. Its numerical orbits exhibited persistent nonperiodic motion and instability to small changes of the initial condition. Lorenz discussed the consequences for very-long-range weather forecasting: a deterministic law need not make a detailed distant future reliably recoverable from approximate observations. The historically relevant lesson of the Lorenz system is this growth of uncertainty under deterministic evolution.^[ref]

For $L=(D_t+a)(D_t+b)$ with $a>b>0$, $H$ denotes the Heaviside step function. This Green function vanishes at negative time and satisfies $LG=\delta (t)$. Its positive-time part solves the homogeneous ordinary differential equation, is zero at the origin, and has right derivative one. Thus it is continuous with a unit derivative jump, which gives precisely the Dirac delta function in its second distributional derivative.

With the negative-exponential Fourier transform convention,Both poles lie in the upper half-plane. Contour integration closes above for positive time and below for negative time, giving the title expression. When the integral converges, convolution $y(t)={\int }_{-\infty }^{t}G(t-s)f(s)ds$ supplies the response with no freely added homogeneous transient growing toward negative infinity. The continuous limit $a\to b$ is $G(t)=H(t)t{e}^{-bt}$.

CDK1 activity and protein abundance are different observables. A total-CDK1 Western blot estimates protein amount under its calibration and loading assumptions. A phosphosite-specific Western blot asks about a particular CDK1 phosphorylation state. A histone H1 kinase assay asks how much substrate modification occurs under specified assay conditions. Any conversion between them needs a model and suitable controls.

An original counting example stipulates $100$ CDK1 molecules, with $60$ bound to cyclin B in each of two hypothetical samples. In sample A, $10$ bound molecules are activity-permissive and $50$ are inhibited; in sample B, $40$ are activity-permissive and $20$ are inhibited. The other $40$ are unbound in both samples. Total CDK1 and bound CDK1 agree, but activity-permissive counts differ fourfold. These are invented counts for bookkeeping, not measured biological populations.

If each activity-permissive molecule contributes the same reaction rate, substrate and ATP remain available, and the assay stays linear in time, a simple illustrative readout isHere $R$ is the recorded signal, $R_0$ is a matched blank signal, $N_{\mathrm{active}}$ is the activity-permissive count, $t$ the reaction interval, and $\alpha$ a common signal contribution per active molecule per unit time. The corrected signals in the stipulated samples have ratio $40/10=4$, although their total protein amounts agree. This is an assay accounting assumption, not an ordinary differential equation for the cell cycle. The no-added-cyclin B baseline in the inspected experiment can itself contain catalytic activity; it need not equal an enzyme-free blank $R_0$.

Substrate depletion, differing inhibitory conditions or a saturated detector can break the simple proportionality. Comparing samples therefore requires matched reaction intervals, substrate supply, loading and detector response, rather than assuming that brighter histone H1 labeling automatically means more total CDK1.

The characteristic ordinary differential equations $\dot{x}=y$, $\dot{y}=x$ preserve $u=x^2-y^2$. Away from the origin, $u=x^2-y^2$, $v=2xy$ have nonzero Jacobian determinant $4(x^2+y^2)$, and the chain rule reduces the displayed equation to $2(x^2+y^2)f_v=0$. Hence locally $f=F(u)$.

Globally, a level set of $u$ can have disconnected characteristic components, permitting different profiles on its different branches. A smooth counterexample on the whole plane is $f=\exp [-1/(x^2-y^2{)}^2]$ where $x>0$ and $x^2>y^2$, and $f=0$ elsewhere. The exponential is flat on the wedge boundary, so this extension is smooth and satisfies the equation. Nevertheless the points $(1,0)$ and $(-1,0)$ have equal $u$ and different function values. One global profile $F(u)$ therefore needs an additional branch-identification condition, rather than following from the local change of variables alone.

For $c(y)=1-y^2$ in the strip $|y|<1$, the travel-time coordinate is $s(y)=\mathrm{artanh}y$. The characteristic curves areequivalently $y=\tanh (C\mp x)$. Outside the strip the same separated ordinary differential equation gives hyperbolic cotangent branches, and the equilibrium curves $y=\pm 1$ are also characteristic.

For the first-order quasilinear partial differential equation $P(x,y,z)z_x+Q(x,y,z)z_y=R(x,y,z)$, the displayed ordinary differential equations trace curves tangent to every classical solution surface. Indeed, if $z=u(x,y)$ solves the equation and $x_s=P$, $y_s=Q$, then the chain rule gives $z_s=u_{x}P+u_{y}Q=R$. Conversely a surface swept out by these curves solves the equation wherever its projection to $(x,y)$ is invertible and hence represents a differentiable graph. The coefficients can depend on the unknown $z$, which distinguishes these equations from characteristic equations for a prescribed linear transport equation.

A characteristic speed is an eigenvalue of the flux Jacobian. Along a corresponding characteristic curve, an associated wave variable obeys an ordinary differential equation obtained by projecting the governing system onto the matching left eigenvector.

At chemical equilibrium, a reversible chemical reaction has equal forward and reverse reaction rates, so its macroscopic composition stays constant. Individual molecules can keep changing identity: the reaction has not stopped. This dynamic balance is the definition described by the ^[ref]. Equal reaction rates also need not mean equal concentrations. A small population converting rapidly can balance a large population converting slowly.

This distinction gives two useful questions. Which composition is favored at the specified temperature and other constraints? Thermodynamics answers through chemical potentials and an equilibrium constant. How quickly does a prepared mixture approach that composition? Chemical kinetics answers through reaction rates. Knowing the destination does not determine the journey's duration. A composition that barely changes during observation may reflect extremely slow chemical kinetics, so constancy alone does not establish chemical equilibrium.

Consider the idealized elementary chemical reactionAssume a well-mixed closed system at fixed temperature and volume, with first-order law of mass action rates in both directions and no competing chemical reactions. Write $a=[\mathrm{A}]$ and $b=[\mathrm{B}]$ for concentrations in $\mathrm{M}=\mathrm{mol}/\mathrm{L}$. The total $a+b=C=1\mathrm{M}$ is conserved because each conversion replaces one molecule by one molecule. This conservation is special to the chosen stoichiometry.

The forward and reverse reaction rates, measured as conversions per unit volume per unit time, are $J_f=k_{f}a$ and $J_r=k_{r}b$. HenceSetting this derivative to zero givesWith $a(0)=1\mathrm{M}$ and $b(0)=0$, solving the ordinary differential equation yieldsAt $t=5\mathrm{s}$, these are approximately $0.3515\mathrm{M}$ and $0.6485\mathrm{M}$. Initially $J_f=0.3M/s$ and $J_r=0$. At chemical equilibrium, both are $0.075M/s$, while the net reaction rate $J_f-J_r$ is zero. The persistent exchange is visible in the lower panel even after the upper panel has nearly flattened.

Why does chemical equilibrium also involve a free energy minimum? For a closed system at fixed temperature and pressure, without imposed non-expansion work, stable equilibrium minimizes Gibbs free energy over compositions allowed by conserved elemental inventories. See W. C. Carter, MIT thermodynamics notes, sections 2-1.14 and 2-1.19. This is constrained Gibbs minimization for chemical equilibrium: freely minimizing over unrelated species amounts would change the inventory. At fixed temperature and volume, the appropriate potential is Helmholtz free energy. In either setting, varying an allowed extent of reaction gives the corresponding free energy slope ${\sum }_i{\nu }_i{\mu }_i$, with molar chemical potentials and signed stoichiometric vector entries. An interior equilibrium has zero slope; stability requires a minimum. For A converting to B, this says ${\mu }_{\mathrm{B}}={\mu }_{\mathrm{A}}$.

The beauty of this connection is that molecular exchange, a macroscopic composition, and a geometric minimum describe the same balance. In the ideal concentration model, let $x=b/C$ and choose the standard molar chemical potential difference consistent with $K=3$. Up to a composition-independent constant, the fixed-volume Helmholtz free energy then obeyswhere $n$ is the conserved total amount of substance and $R$ the molar gas constant. Its derivative is $\ln [x/(1-x)]-\ln 3$, vanishing at $x=3/4$. Its second derivative is $1/x+1/(1-x)>0$ inside $0<x<1$, so this balance is the unique minimum. The mixing entropy terms explain why the favored composition still contains A.

A chemostat maintains a well-mixed culture at constant volume $V$ by feeding fresh nutrient solution and removing culture at the same volumetric flow $Q$. Its dilution rate is $D=Q/V$. Mixing makes the outflow's nutrient concentration and biomass concentration equal to those inside the vessel; continuous replacement supplies nutrient and also removes growing cells.

The simplest one-nutrient chemostat uses nutrient concentration $s(t)$, biomass concentration $x(t)$, constant feed nutrient concentration $s_{\mathrm{in}}>0$, and constant biomass yield $Y>0$. Assume a cell-free feed, one growth-limiting nutrient, constant volume, and no cell retention, death or additional nutrient consumption for maintenance. With specific growth rate $\mu (s)$, its ordinary differential equations areThe nutrient balance contains inflow $D{s}_{\mathrm{in}}$, outflow $Ds$, and consumption $\mu (s)x/Y$; the biomass balance contains production $\mu (s)x$ and removal $Dx$. If $s$ is in grams of nutrient per litre and $x$ in grams of dry biomass per litre, $Y$ is grams of dry biomass produced per gram of nutrient consumed. Thus both balances have units of their respective concentration per time. The Monod equation supplies an empirical choice of $\mu (s)$; the balances themselves express the stated production, consumption and transport assumptions.

The appeal of a chemostat is that, at a positive chemostat steady state, adjusting the pump sets the specific growth rate to $D$. The residual nutrient concentration then adjusts to that rate, and the available nutrient together with biomass yield sets the biomass concentration. This gives a controlled way to compare cultures at selected growth rates and a simple quantitative warning of chemostat washout.

Two 1950 publications on continuous culture are Novick and Szilard's Description of the Chemostat, Science 112(2920), 715-716, dated 15 December, and Monod's La technique de culture continue. Théorie et applications, Annales de l'Institut Pasteur 79, 390-410. The latter's original bibliographic details are reproduced in the contents of Monod's collected papers. These references identify contemporary publications without establishing priority between them. ^[ref]. ^[ref].