Singular Perturbation Derivation of a Reduced ODE

Why a Small Parameter Produces a Reduced ODE

A QSP model can contain states that change on substantially different timescales. When one state responds much faster than the others, it may rapidly approach a relation determined by the slower states. The resulting approximation can replace the full system with an ODE for the slow states alone. The key issue is that this replacement is justified only when the fast dynamics actually approach the selected relation, including from the specified initial condition.

Write a system with slow state \(x\in\mathbb{R}^m\), fast state \(y\in\mathbb{R}^q\), and dimensionless small parameter \(0<\varepsilon\ll1\) as

\frac{dx}{dt}=f(x,y,\varepsilon),\qquad \varepsilon\frac{dy}{dt}=g(x,y,\varepsilon).

Here \(t\) is slow time. The factor \(\varepsilon\) makes the \(y\) dynamics fast relative to \(x\), provided the relevant rates in \(f\) and \(g\) are of comparable order after scaling. This form is typically obtained after nondimensionalization; a small numerical coefficient alone is not enough if the state and time scales have not been interpreted consistently.

Definition: A leading-order singular-perturbation reduction is obtained by finding a branch \(y=h(x)\) satisfying \(g(x,h(x),0)=0\), then using the reduced slow ODE \(dx_r/dt=f(x_r,h(x_r),0)\). The branch is an approximation to the fast state only where the fast dynamics attract trajectories to it.

Deriving the Leading-Order System

First set \(\varepsilon=0\) in the fast equation. This does not mean setting \(dy/dt=0\) in the original equation as an exact identity. Instead, it identifies the constraint that the fast state must satisfy to remain near the slow evolution at leading order:

$$ g(x,y,0)=0. $$

Suppose this equation defines a sufficiently smooth, single-valued branch \(y=h(x)\) in the region of interest. Substituting that branch into the slow equation gives

\frac{dx_r}{dt}=f(x_r,h(x_r),0),\qquad x_r(0)=x_0.

The subscript \(r\) denotes the reduced solution. The usual leading-order initial condition is the original slow initial state \(x_0\); the initial fast state \(y_0\) is not forced to equal \(h(x_0)\). That mismatch is handled by a short initial layer, rather than by changing the leading-order slow initial condition.

A formal expansion clarifies the approximation. Away from the initial layer, suppose \(y=h(x)+\varepsilon q+O(\varepsilon^2)\). Write \(g_y\) and \(g_\varepsilon\) for derivatives evaluated at \((x,h(x),0)\), and let \(f_0=f(x,h(x),0)\). At order \(\varepsilon\), the fast equation gives

D h(x)f_0=g_yq+g_\varepsilon,\qquad q=g_y^{-1}\bigl(D h(x)f_0-g_\varepsilon\bigr),

when \(g_y\) is invertible. The leading-order equation is therefore the constraint \(g(x,h(x),0)=0\); the correction \(q\) describes a possible first-order displacement from it. This expansion is an outer approximation and does not describe an initial mismatch of order one.

Initial-Layer Analysis and Conditions for Validity

Introduce fast time \(\tau=t/\varepsilon\). The equations become

\frac{dx}{d\tau}=\varepsilon f(x,y,\varepsilon),\qquad \frac{dy}{d\tau}=g(x,y,\varepsilon).

On a bounded interval of \(\tau\), \(x\) changes only by order \(\varepsilon\), so to leading order the initial-layer trajectory \(Y(\tau)\) solves

\frac{dY}{d\tau}=g(x_0,Y,0),\qquad Y(0)=y_0.

For the proposed reduction to apply, this trajectory must approach \(h(x_0)\). Local stability of the branch by itself does not establish that result: a locally stable equilibrium attracts states in its basin, not every possible initial fast state. As \(x\) evolves, attraction must also persist for the slow states encountered over the time interval.

Theorem: A standard finite-time approximation conclusion holds under regularity and boundedness conditions such as the following: \(f\) and \(g\) are sufficiently smooth on a compact region containing the relevant solutions; \(g(x,h(x),0)=0\) defines a smooth branch there; the eigenvalues of \(g_y(x,h(x),0)\) have real parts uniformly bounded away from zero on the stable side; and solutions remain in that region. In addition, the initial fast state must lie in a suitable basin of attraction of \(h(x_0)\), and its initial-layer trajectory must be attracted to that branch uniformly over the relevant slow states. With these conditions, for fixed \(T\) on which the solutions remain in the region, the slow-state error is typically bounded by \(C_T\varepsilon\) on \(0\le t\le T\), with \(C_T\) independent of sufficiently small \(\varepsilon\). The fast-state error need not be small at \(t=0\); after the initial layer it is typically order \(\varepsilon\) under the same smoothness and stability conditions.

The slow-state error can be order \(\varepsilon\) despite an order-one initial mismatch in \(y\): the mismatch acts over a time of order \(\varepsilon\), so its accumulated effect on \(x\) can be small. This reasoning depends on the initial layer being attracted to the chosen branch and on its effect being integrable and bounded. A solution that escapes the branch’s basin does not satisfy this argument.

Worked Example: Exact Initial-Layer Error in a Scalar System

Consider the dimensionless system \(x'=y-x\), \(\varepsilon y'=x-y\), with \(x(0)=1\) and \(y(0)=0\). Here a prime denotes differentiation with respect to \(t\). The leading constraint is \(x-y=0\), so \(h(x)=x\). Substitution gives the reduced ODE \(x_r'=x_r-x_r=0\), with solution \(x_r(t)=1\).

To check the approximation, define \(d=y-x\). The full system gives \(x'=d\) and \(y'=-d/\varepsilon\), so

d'=y'-x'=-\left(\frac{1}{\varepsilon}+1\right)d,\qquad d(0)=-1.

Thus \(d(t)=-\exp[-(1+1/\varepsilon)t]\). Integrating \(x'=d\) yields

x(t)=1-\frac{\varepsilon}{1+\varepsilon} \left(1-\exp\left[-\frac{1+\varepsilon}{\varepsilon}t\right]\right).

Consequently, \(|x(t)-x_r(t)|\le\varepsilon/(1+\varepsilon)\) for all \(t\ge0\), while \(y(0)-h(x(0))=-1\). The fast discrepancy decays on a time scale of order \(\varepsilon\), and the slow-state error remains order \(\varepsilon\). The fast equilibrium \(y=x\) is attracting for every fixed \(x\), and the initial \(y=0\) lies in its basin.

A Stable Branch Is Not Enough

The next example shows why the attraction condition on the initial fast state is essential. Consider \(x'=y\) and \(\varepsilon y'=y-y^3\). The fast equation has equilibria \(y=1\), \(y=0\), and \(y=-1\). Since the derivative of \(g(y)=y-y^3\) is \(1-3y^2\), the equilibria \(y=1\) and \(y=-1\) are stable, while \(y=0\) is unstable.

Worked Example: Initial State Outside the Selected Branch’s Basin

Select the stable branch \(h(x)=1\), and set \(x(0)=0\), \(y(0)=-1\). The formal reduced equation is \(x_r'=1\), giving \(x_r(t)=t\). But \(y=-1\) is itself an equilibrium of the full fast equation, so the full solution has \(y(t)=-1\) and \(x(t)=-t\). Therefore

$$ |x(t)-x_r(t)|=2t. $$

For any fixed positive \(t\), this discrepancy does not tend to zero as \(\varepsilon\) decreases. The branch \(h=1\) is locally stable, but the initial state \(y_0=-1\) lies in the basin of the other stable branch. The reduction based on \(h=1\) is therefore invalid for these initial conditions. A valid analysis must select a branch that attracts the initial-layer trajectory, or restrict the initial conditions to the selected branch’s basin.

Worked Reduction in a Two-State Slow System

Worked Example: Propagating the Error Through a Slow State

Let \(x_1\) and \(x_2\) be dimensionless slow states and \(y\) a dimensionless fast state:

x_1'=-x_1+y,\qquad x_2'=x_1,\qquad \varepsilon y'=x_1-y.

For each fixed \(x_1\), the fast constraint gives \(h(x_1,x_2)=x_1\); its fast derivative is \(-1\), so the branch is uniformly attracting. Take \((x_1(0),x_2(0),y(0))=(2,1,0)\), which is in the branch’s basin. The reduced equations and solution are

x_{1,r}'=0,\qquad x_{2,r}'=x_{1,r}, \qquad x_{1,r}(t)=2,\qquad x_{2,r}(t)=1+2t.

For the full system, set \(d=y-x_1\). Then \(x_1'=d\), \(y'=-d/\varepsilon\), and \(d'=-(1+1/\varepsilon)d\), with \(d(0)=-2\). Writing \(k=(1+\varepsilon)/\varepsilon\), integration gives

x_1(t)=2-\frac{2\varepsilon}{1+\varepsilon}(1-e^{-kt}),\qquad x_2(t)=1+2t-\frac{2\varepsilon}{1+\varepsilon} \left(t-\frac{1-e^{-kt}}{k}\right).

At \(\varepsilon=0.02\) and \(t=0.5\), \(k=51\) and \(e^{-25.5}\) is negligible at the shown precision. Thus \(x_1(0.5)\approx1.9608\), compared with the reduced value \(2\). For the second state,

x_2(0.5)\approx 2- \frac{0.04}{1.02}\left(0.5-\frac{1}{51}\right) \approx1.9812.

The reduced value is \(2\), so the error in \(x_2\) is about \(0.0188\). This is small but not zero: the brief initial discrepancy in \(x_1\) affects the subsequent evolution of \(x_2\). The calculation also illustrates why checking only the fast-state timescale, without checking how the fast state enters the slow equations, can understate the reduction error.

Diagnostics, Alternatives, and Limitations

A reduction should be assessed over the state and parameter range in which it will be used, not just at one nominal trajectory. Useful checks include:

Fast stability along the branch: inspect the eigenvalues of \(g_y(x,h(x),0)\) across the relevant slow-state range. A loss of uniform attraction can invalidate the approximation.
Initial-layer attraction: simulate or analyze \(dY/d\tau=g(x_0,Y,0)\) from the actual \(y_0\). Confirm that it approaches the selected branch rather than another equilibrium, a cycle, or an unbounded trajectory.
Branch and solution validity: check that the constraint defines the intended smooth branch and that both full and reduced trajectories remain in the region where the assumptions hold.
Full-model comparison: at representative small values of \(\varepsilon\), compare slow-state trajectories with the reduced solution. A persistent order-one discrepancy, abrupt branch departure, or error that fails to decrease as \(\varepsilon\) decreases signals trouble.

Directly integrating the full system retains the fast transient and is the main alternative when that transient matters or the assumptions are uncertain, though the separated timescales may make numerical integration demanding. The reduced ODE is cheaper and easier to analyze, but it omits the initial fast transient and its finite contribution to slow-state changes. Simply substituting an equilibrium relation without verifying branch regularity, stability, and basin attraction is not a sound alternative; it is an unvalidated shortcut.

The finite-time error statement is not a guarantee over arbitrarily long times. Small errors can accumulate, and parameter changes or evolving states can move the system toward a region where attraction weakens or the branch ceases to be appropriate. The leading-order reduction also does not preserve all details of the fast trajectory. If those details are important to the scientific question, the full system or a higher-order approximation may be needed.

Common mistake: “The selected fast equilibrium is stable, so any initial fast state approaches it and the reduced ODE is valid.” Stability is generally local; the initial state may lie outside that equilibrium’s basin. Instead, verify that the initial-layer trajectory from the specified \(y_0\) is attracted to the selected branch, and that attraction persists over the slow states of interest.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

For \(x'=f(x,y,\varepsilon)\), \(\varepsilon y'=g(x,y,\varepsilon)\), what equation defines the leading-order fast constraint?
Why does the leading reduced ODE usually use \(x_r(0)=x_0\) even when \(y_0\ne h(x_0)\)?
A branch is locally stable, but the initial fast state belongs to the basin of a different stable equilibrium. Does the usual finite-time approximation conclusion follow?
In the first worked example, what is the maximum slow-state error for \(t\ge0\), and how does it scale with \(\varepsilon\)?
Name two checks that can reveal a failure of the reduction before it is used for interpretation.

Solutions. (1) Solve \(g(x,h(x),0)=0\) for the selected branch \(h(x)\), then use \(x_r'=f(x_r,h(x_r),0)\). (2) The slow state changes by order \(\varepsilon\) during a bounded fast-time initial layer; the fast mismatch is treated through that layer, not by replacing the leading slow initial condition. (3) No. Stability of the branch does not establish attraction from an initial state outside its basin. (4) The error is bounded by \(\varepsilon/(1+\varepsilon)\), which is order \(\varepsilon\) for small \(\varepsilon\). (5) Examples include checking fast-subsystem attraction from the actual initial state, examining fast Jacobian eigenvalues along the branch, verifying that the branch and trajectories remain in the valid region, and comparing full and reduced trajectories as \(\varepsilon\) decreases.