Numeric Integration for DAEs

Prepared by: Prof. Alexander Dowling, Myia Dickens (mdicken2@nd.edu,2023)

This notebook is the computational companion to Lectures 5--7. It connects error and stability experiments for ODE integrators to index-1 DAE integration, then prepares for direct collocation in Pyomo.dae.

Reference guide¶

The lecture develops the main derivations through partial fill-ins. This guide retains the detailed assumptions, calculations, and implementation choices for later use.

• Quadrature, precision, and accumulated error

• Why Gauss nodes give precision $2L$

• From an integral to a time step

• Local error, global error, and zero-stability

• Higher-order explicit methods

• Scalar stability and oscillation

• Complex stability regions and stiff decay

• Time-varying contraction and algebraic stability

• Stiffness and solver families

• Step changes, control jumps, and method selection

• RK stages, tableaus, and stability functions

• Index-1 DAE stages and coupled residuals

• Integration within dynamic optimization

• Sources and further reading

Quadrature, precision, and accumulated error¶

For $h=b-a>0$, approximate an integral by function samples:

The nodes $c_j$ lie in $[0,1]$ for the rules considered here. The weights are dimensionless. A rule has precision $p$ if it integrates every polynomial of degree at most $p-1$ exactly, and fails for at least one degree-$p$ polynomial. This convention makes precision one larger than the maximum degree of exactness.

The bounds require the relevant derivatives to remain bounded as the mesh is refined. For example, expand about the midpoint $m=(a+b)/2$:

Odd powers cancel by symmetry. Midpoint quadrature omits the second term. Trapezoidal quadrature instead gives

so $I-Q_{\rm trap}=-h^3 f''(m)/12+O(h^5)$. These are single-interval statements. On a fixed interval of length $T$, a uniform composite rule sums $N=T/h$ interval errors. If each error is bounded by $C h^{p+1}$ with the same $C$, then

This is why the composite rule loses one order. For nonuniform widths $h_i\le h_{\max}$, the same argument gives $\sum_i C h_i^{p+1}\le CT h_{\max}^{p}$. Cancellation can improve a particular example; the bound does not assume it.

The forward finite difference is $[f(b)-f(a)]/(b-a)$. Reversing the numerator changes its sign. Quadrature accuracy and finite-difference accuracy concern different approximations and should not be interchanged. Ascher & Petzold, p. 76, supplies the one-interval precision convention; the composite bound above follows by summation.

Why Gauss nodes give precision $2L$¶

One node on the symmetric interval¶

Require $Q(f)=w f(c)$ to integrate every linear function on $[-1,1]$ exactly. Matching the constant and linear moments gives

Thus $c=0$ and $Q(ax+b)=2b=\int_{-1}^{1}(ax+b)\,dx$. The rectangle at the sampled height has equal excess and deficit areas: the odd contribution $ax$ cancels. The sketch uses $a=1/2,b=1$; each triangle has area $1/4$. Reversing the slope swaps their roles. The rule fails for $f(x)=x^2$, since $Q(f)=0$ while the integral is $2/3$, so its precision is two.

Orthogonal polynomial construction on the unit interval¶

Let $P_L$ be a degree-$L$ polynomial orthogonal to every lower-degree polynomial on $[0,1]$ with unit weight:

Choose its $L$ distinct roots $c_1,\ldots,c_L$ as nodes. With the Lagrange basis

the quadrature is exact on $\mathcal P_{L-1}$. Now take any $p\in\mathcal P_{2L-1}$. Polynomial division gives

At every node, $p(c_j)=r(c_j)$ because $P_L(c_j)=0$. Orthogonality also gives $\int_0^1qP_L=0$. Therefore

Both zero terms matter: roots remove $qP_L$ from the discrete sum, while orthogonality removes it from the integral. Counting $L$ node locations and $L$ weights suggests $2L$ matching moments but does not prove that suitable real nodes exist or that all moment equations can be satisfied. The orthogonal-polynomial construction does.

No real $L$-node formula can be exact through degree $2L$: the nonnegative polynomial $\prod_j(\tau-c_j)^2$ vanishes at all nodes but has a positive integral. Thus the maximum precision is $2L$. With $L=1$, $P_1(\tau)=\tau-1/2$, so $c_1=1/2$ and $w_1=1$. With $L=2$, the nodes are $1/2\pm1/(2\sqrt3)$ and the weights are $1/2,1/2$.

Biegler, Theorem 10.1, p. 291, uses this construction. Here “degree at most $2L-1$” states the convention without relying on the book’s ambiguous use of polynomial “order.”

From an integral to a time step¶

For $\dot{\mathbf z}=\mathbf f(t,\mathbf z)$, the exact identity is

The integrand contains the unknown trajectory. Using the left endpoint gives forward Euler (FE), the right endpoint gives backward Euler (BE), and averaging both gives the trapezoidal method, also called Crank-Nicolson (CN):

FE evaluates a known state. BE and CN require a solve for the new state. A smaller quadrature error does not remove the cost or possible failure of that nonlinear solve. Ascher & Petzold, pp. 77-78, derives these formulas.

Local error, global error, and zero-stability¶

Write a difference operator $\mathcal N_h$ with the derivative normalization included. For FE,

The normalized local truncation error is $\mathbf d_i=(\mathcal N_h\mathbf z)_i$, obtained by inserting the exact trajectory. Taylor expansion gives $\mathbf d_i=\tfrac12h_i\mathbf z''(t_i)+O(h_i^2)$. Thus FE is consistent of order one. The corresponding unnormalized one-step defect is $h_i\mathbf d_i=O(h_i^2)$. Some references call that defect the local error; state the normalization whenever quoting an order.

The global error at a mesh point is $\mathbf e_i=\mathbf z(t_i)-\mathbf z_i$, using a fixed initial-value problem and a fixed final time. Convergence of order $q$ means $\max_i\|\mathbf e_i\|\le C h_{\max}^{q}$ for sufficiently small $h_{\max}$, with $C$ independent of the mesh.

Zero-stability controls the propagation of perturbations as the mesh is refined. One operator form, on the admissible nearby grid functions, is

where $K$ and the allowed maximum step are independent of the mesh on a fixed time interval. This is a statement about a family of shrinking meshes, rather than a large-step absolute-stability condition.

To see consistency plus zero-stability imply convergence, substitute the exact grid values for $\mathbf x$ and the computed solution for $\mathbf v$. The computed residual is zero; the exact residual is $\mathbf d$. Hence

If initial error and normalized truncation error are $O(h_{\max}^{q})$, so is the global error. This proof requires the norm estimate, regularity, and fixed time interval; “add $N$ errors” alone does not establish stability for an arbitrary method. Ascher & Petzold, Definition 3.2 and Theorem 3.1, pp. 38-39; Biegler, equations (9.5)-(9.6), p. 254.

Method orders and error normalization¶

Under the smoothness and zero-stability assumptions on a fixed interval:

The order concerns refinement toward zero step size. It does not imply that a coarse step is stable or resolves a fast transient. Ascher & Petzold, pp. 77-80, gives the formulas; the preceding operator estimate connects their normalized local orders to the global orders.

Higher-order explicit methods¶

The explicit midpoint rule estimates the state halfway through a step, then uses its derivative for the full update:

It has normalized local order two and one-step defect $O(h_i^3)$. The classical four-stage RK4 formula is

Its normalized local order is four and its one-step defect is $O(h_i^5)$. These global orders require sufficient smoothness and stability. Higher order reduces the small-step error but does not imply A-stability. On $z'=\lambda z$, explicit midpoint has $R(\zeta)=1+\zeta+\zeta^2/2$, while RK4 has $R(\zeta)=1+\zeta+\zeta^2/2+\zeta^3/6+\zeta^4/24$; both have bounded absolute-stability regions. Ascher & Petzold, equations (4.3) and (4.7), pp. 78-79.

Single-Step Runge-Kutta Methods¶

• Chapter 9 in Biegler (2010)

• Chapter 17 in McClarren (2018)

General Form: Index 0 DAE¶

Consider the ODE system:

where $\mathbf{z}(t)$ are the differential variables and $\mathbf{f}(t,\mathbf{z})$ is a (nonlinear) continuous function.

The general Runge-Kutta formula is:

where

• $\mathbf{z}_i$ are the differential variables at the start of step $i$ (time $t_i$)

• $\mathbf{z}_{i+1}$ are the differential variables at the end of step $i$ (time $t_{i+1}$)

• $\hat{\mathbf{z}}_k$ are differential variables for intermediate stage $k$

• $h_i$ is the size for step $i$ such that $t_{i+1} = t_i + h_i$

• $n_s$ is the number of stages

• $a_{k,j}$ are coefficients, together known as the Runge-Kutta matrix

• $b_k$ are coefficients, known as the weights

• $c_k$ are coefficients, known as the nodes

$h_i$ is selected using local error estimates and stability requirements. Zeros in $\mathbf A$ express which stage evaluations are absent.

The choice for $\mathbf{A}$, $\mathbf{b}$ and $\mathbf{c}$ selects the specific method in the Runge-Kutta family. These coefficients are often specified in a Butcher block (or Butcher tableau).

The first-order condition is $\mathbf b^T\mathbf 1=1$. For the internally consistent stage times used here, $\mathbf A\mathbf1=\mathbf c$ as well:

Explicit (Forward) Euler¶

Consider one of the simplest Runge-Kutta methods:

What are $\mathbf{A}$, $\mathbf{b}$ and $\mathbf{c}$ in the general formula?

$n_s = 1$. This is only a single stage. Thus we only need to determine $c_1$, $b_1$, and $a_{11}$

Moreover, $\hat{\mathbf{z}}_1 = \mathbf{z}_i$ because $\mathbf{f}(\cdot)$ is only evaluated at $t_i$ and $\mathbf{z}_i$. This implies:

• $a_{11} = 0$

• $c_1 = 0$

• $b_1 = 1$

• $\mathbf{A}=[0]$; this is a one-by-one zero matrix, not an empty matrix.

The implementation is very straightforward (see below). We can calculate $\mathbf{z}_{i+1}$ with a single line!

Implicit (Backward) Euler¶

Consider another simple Runge-Kutta method:

What are $\mathbf{A}$, $\mathbf{b}$ and $\mathbf{c}$ to express using the general formula?

$n_s = 1$. This is only a single stage. Moreover, $\hat{\mathbf{z}}_1 = \mathbf{z}_{i+1}$ because $\mathbf{f}(\cdot)$ is evaluated at $t_{i+1}$ and $\mathbf{z}_{i+1}$. This implies:

• $b_1 = 1$

• $c_1 = 1$

Moreover, $\mathbf{z}_{i+1} = \mathbf{z}_{i} + h_i~\mathbf{f}(t_{i+1}, \mathbf{z}_{i+1})$ implies $\hat{\mathbf{z}}_1 = \mathbf{z}_i + h_i \mathbf{f}(t_{i+1}, \hat{\mathbf{z}}_1)$. Thus:

• $a_{1,1} = 1$

Notice that the formula for $\mathbf{z}_{i+1}$ is implicit. We need to solve a (nonlinear) system of equations to calculate the step.

Key Differences¶

Comparison¶

Let’s test this on a simple problem:

The solution to this problem is

For simplicity, let’s numerically analyze $\lambda = 1$.

The cells below compute both panels. Their common plotting cell is also the source of the handout image. Change the step sizes and rerun to test another mesh.

Stability¶

Keeping $\lambda = 1$, are there any limits on step size?

Scalar stability and oscillation¶

For $\dot z=-\alpha z$, $\alpha>0$, $z_0>0$, set $\mu=\alpha h$. The exact solution is positive and decreasing. Each numerical method gives $z_{i+1}=g(\mu)z_i$. Substitute the test equation into each update and solve for the new state:

For this scalar problem, $|g|\le1$ means no amplification. The stronger condition $0\le g\le1$ also preserves nonnegative monotone decay. This sign test is a direct observation about the scalar recurrence, not a general nonlinear stability theorem.

At FE’s endpoint $\mu=2$, values alternate with unchanged magnitude. They are bounded, but do not approach zero. CN can also alternate when $\mu>2$ even though its magnitude decreases. Thus stability, non-oscillation, and accurate transient resolution are separate questions. The prepared comparison uses FE and BE at $h=1$ and $h=2.5$ with $\alpha=1$.

Complex stability regions and stiff decay¶

Apply a method to $\dot z=\lambda z$ with complex $\lambda$ and let $\zeta=h\lambda$. If $z_{i+1}=R(\zeta)z_i$, its absolute-stability region is the set where $R$ is defined and $|R(\zeta)|\le1$. For FE, BE, and CN:

FE’s region is the disk $|1+\zeta|\le1$. BE’s is the exterior $|1-\zeta|\ge1$. For CN, squaring $|1+\zeta/2|\le|1-\zeta/2|$ gives $\operatorname{Re}\zeta\le0$. A method is A-stable if its region contains the entire closed left half-plane, with no poles there. BE and CN satisfy this definition; FE does not.

No explicit RK method can be A-stable: its nonconstant polynomial $R(\zeta)$ is unbounded as $\zeta\to-\infty$. For a linear vector system, eigenvalues provide scalar mode tests when such a decomposition is applicable; nonnormal coupling can cause transient growth that scalar eigenvalues alone do not describe. If $\operatorname{Re}\lambda=0$ and $\operatorname{Im}\lambda\ne0$, the exact scalar solution oscillates at constant magnitude. When $\lambda=0$, it is constant.

For a bounded smooth moving equilibrium $\gamma(t)$, consider

The stiff-decay requirement used by Biegler is that the computed endpoint follows the moving equilibrium, $|z_i-\gamma(t_i)|\to0$ as $\operatorname{Re}(h\lambda)\to-\infty$, at fixed positive step times after the initial layer. In the homogeneous special case $\gamma=0$, this requires $R(\zeta)\to0$. BE passes that test and CN fails: CN’s limit is -1, so a very fast decaying mode can remain as alternating numerical values.

The homogeneous limit alone is not a complete test for an arbitrary forced RK formula. A sufficient structural condition for the implicit RK setting discussed by Ascher & Petzold is nonsingular $\mathbf A$ together with stiff accuracy ($\mathbf b^T$ equals the final row of $\mathbf A$, with endpoint stage $c_s=1$). It makes the last stage the new state, tying the endpoint to the equilibrium. Standard Radau IIA methods have this property. Standard Gauss methods and the Lobatto IIIA family used for trapezoidal collocation do not have stiff decay. Ascher & Petzold, pp. 42-43, 56-57, 102-104; Biegler, pp. 255-256.

Time-varying contraction and algebraic stability¶

A-stability tests a constant $\lambda$. AN-stability tests $z'=\lambda(t)z$ with $\operatorname{Re}\lambda(t)\le0$. For an $s$-stage RK method, let $\mathbf Z=\operatorname{diag}(\zeta_1,\ldots,\zeta_s)$, where $\zeta_j=h\lambda(t_i+c_jh)$. The scalar amplification is

AN-stability requires this function to be defined and to satisfy $|K(\mathbf Z)|\le1$ for every admissible stage tuple with $\operatorname{Re}\zeta_j\le0$. Repeated nodes must use equal $\zeta_j$, because they sample the same time. When all $\zeta_j=\zeta$, $K$ reduces to $R(\zeta)$.

The real-coefficient method is algebraically stable if

Equivalently, $M_{ij}=b_i a_{ij}+b_j a_{ji}-b_i b_j$. Algebraic stability implies AN-stability. The converse holds when the nodes $c_j$ are distinct. Burrage & Butcher (1979), Theorem 3.3, printed p. 50, states that qualification explicitly; the unqualified textbook synonym is too broad for repeated-node RK formulas.

For BE, $A=b=c=1$, so $M=1\ge0$. For FE, $A=c=0$, $b=1$, so $M=-1$, which fails. For trapezoidal RK, $\mathbf A=[[0,0],[1/2,1/2]]$, $\mathbf b=(1/2,1/2)^T$, $\mathbf c=(0,1)^T$, giving $\mathbf M=\operatorname{diag}(-1/4,1/4)$: A-stability does not imply algebraic stability. Gauss and Radau IIA collocation have algebraic stability; this does not give Gauss stiff decay. These properties answer different questions.

For a direct time-varying counterexample, trapezoidal RK has

Choose $\zeta_1=-10$ and $\zeta_2=0$, both in the left half-plane. Then $K=-4$, so the stage-varying problem can amplify even though the constant-coefficient function $R(\zeta)=K(\zeta,\zeta)$ is A-stable. This separates the two tests without relying on an unexplained alias.

Key observation: forward (explicit) Euler becomes unstable with large steps whereas backward (implicit) Euler is stable.

The Lecture 6 handout derives the growth factors and then generalizes them to absolute-stability regions.

Key results (for this specific test problem):

• Forward Euler is stable for this problem when $h\lambda \leq 2$ and non-oscillatory when $h\lambda \leq 1$.

• Backward Euler is stable and non-oscillatory for every $h > 0$ when $\lambda > 0$.

• Crank--Nicolson is stable for every $h > 0$, but oscillates when $h\lambda > 2$.

• Similar analysis and concepts extend to Runge-Kutta methods.

Why the step-size limit exists, in one picture. Rendered from figures/plots/stability-regions.py; the same image printed in the course pack.

The shaded set is where the method does not amplify. A method is A-stable when that set covers the entire left half-plane (Ascher & Petzold, p. 56). Forward Euler’s region is bounded, so it never can be -- hence $h\lambda\le 2$ (boundedness, including the neutral endpoint). Backward Euler and Crank-Nicolson both are.

Error Analysis¶

How does our choice in step size $h$ impact the error of these numerical techniques?

The next calculation measures the RMS error and prints the slope over the two finest meshes. Its plotting cell produces the same image used in the handout. The triangle is a supplied scale reference; estimate each data curve’s slope independently.

Both Euler methods are first order in $\Delta t$: halving $\Delta t$ approximately halves the global error. Here we measure the error with the root-mean-square norm

where $N$ is the number of steps the ODE is solved over.

For both Euler methods, the unnormalized one-step defect is $O(h^2)$, the local truncation error divided by $h$ is $O(h)$, and the global error is $O(h)$ under the smoothness and stability assumptions below. Naming the normalization prevents the same calculation from being called both first and second order.

Index-1 DAE stages and coupled residuals¶

A small genuine index-1 example separates the coupled integration solve from the ODE-constrained optimization example below:

Here $g_y=1$, so the algebraic variable is locally determined by $z$. The exact solution is $z(t)=1/(1+t)$ and $y(t)=1/(1+t)^2$. Backward Euler solves both unknown endpoint values together:

The upper-right sign is positive here because $f=-y$ and the residual subtracts $hf$. Elimination gives $h_i z_{i+1}^2+z_{i+1}-z_i=0$. The root continuous from the positive initial state is

This rearrangement avoids subtracting nearly equal numbers for small $h_i$. With $z_i=1$ and $h_i=1/2$, $z_{i+1}=\sqrt3-1\approx0.73205$, $y_{i+1}\approx0.53590$. The exact differential state at $t=1/2$ is $2/3$; satisfying the algebraic equation exactly does not remove time-discretization error.

For a general regular index-1 DAE, the RK stage equations below enforce $\mathbf g=\mathbf0$ at the stages. They do not automatically enforce an endpoint constraint for a weighted differential endpoint that is not itself a stage. A stiffly accurate method uses its last stage as the endpoint; another method may need an explicit endpoint recovery of the algebraic variables. Consistent initial data and a successful coupled nonlinear solve remain necessary.

Extending Numeric Integration to Index-1 DAEs¶

Consider semi-explicit DAEs:

For a regular semi-explicit index-1 system, $\mathbf g_{\mathbf y}$ is nonsingular along the solution. Enforce the algebraic equations at every stage:

Backward differentiation formulas (BDFs) also extend to index-1 DAEs. Backward Euler is the one-step, first-order member of this family.

At each BDF step, solve the coupled equations for $(\mathbf{z}_{i+1},\mathbf{y}_{i+1})$ with Newton’s method. Using the conventional row-by-variable Jacobian and the residual $\mathbf{z}_{i+1}-\sum_j \alpha_j \mathbf{z}_{i-j+1}-\beta_0 h \mathbf{f}(\mathbf{z}_{i+1},\mathbf{y}_{i+1})=0$, the Jacobian is:

This follows Biegler (2010), equations (9.18)--(9.19), pp. 259--260, with two corrections made explicit: derivatives here use the course’s row-by-variable convention, and the upper-right block has a minus sign when the first residual is written in the form above. The printed equation (9.19) combines the opposite sign in that block with $\mathbf{I}-h\beta_0 \mathbf{f}_{\mathbf{z}}$; direct differentiation shows those signs cannot belong to the same residual, and the published errata do not list the issue.

Order requires a theorem for the method and variable. Regular index-1 structure permits local elimination of $\mathbf y$, but it does not prove unconditional retention of every ODE order or stability property. Gear & Petzold (1984), Theorem 2.1, p. 719, proves fixed-step BDF order $k<7$ with correct initial values and its stated regularity assumptions. This is an order result; it is not a general stability theorem. RK algebraic variables can have a different order from differential variables.

Discussion: Why are implicit methods usually preferred for index-1 DAE systems?

Integration within dynamic optimization¶

The DAE background reference gives the full integration-versus-transcription architecture. In an embedded integration approach, a candidate control or parameter vector is supplied to an IVP solver; the NLP sees a reduced objective and constraints evaluated from the resulting trajectory. In simultaneous transcription, time-indexed states and discretized dynamic residuals are NLP variables and constraints. A large sparse NLP can then exploit local time coupling.

The tradeoff is specific. Embedded integration provides adaptive step control and a dynamically feasible trajectory when integration succeeds, but a poor parameter trial can make integration fail. Repeated sensitivities and path-constraint checks can be costly. Simultaneous transcription exposes sparse derivatives, state bounds, and path constraints at its chosen points; it requires a mesh, suitable initialization, and a check between points. Satisfying collocation constraints alone does not certify the continuous-time path.

The complete ODE-constrained optimization builder below remains useful: it minimizes an accumulated cost $x_3(t_f)$ by choosing $u(t)$, uses three ODEs and a path inequality, and exposes the backward-Euler discretization to Ipopt. It is not an index-1 DAE example. Compare the independent coupled residual example above before attributing algebraic-constraint behavior to this model.

Backward Euler transcription of an ODE-constrained optimization problem¶

Reference: Bynum, Hackebeil, Hart, Laird, Nicholson, Siirola, Watson & Woodruff, Pyomo — Optimization Modeling in Python, 3rd ed., Springer (2021), §12.1, p. 181, equations (12.1)–(12.6). (Earlier versions of this notebook cited “Hart, 2010”, which is the 1st edition — a different book with different pagination.)

The example below is an ODE model with a path inequality, not a DAE. It shows how backward Euler exposes time-indexed states, controls, and constraints to an optimization solver.

Consider the following optimal control problem

Discussion: What variable is the problem attempting to minimize? What variable is being optimized? What types of equations are there?

That t is a bare Python float -- the raw ContinuousSet index that Pyomo passed into the rule. It is not a Pyomo object and it carries no units, so there is nothing to attach units to. This is the cleanest illustration in the course of Pyomo issue #1790 (add units support to ContinuousSet/Pyomo.dae, still open): a ContinuousSet cannot hold units, which is also why a units-correct Pyomo.dae model passes assert_units_consistent before discretization and fails after it.

When a model is physical and needs a time-dependent reference trajectory, the workaround is to have the helper return a bare number and multiply by the intended units at the call site. Hot Air Balloon Dynamic Control does exactly that.

Stiffness and solver families¶

A practical symptom of stiffness is that a stable computation needs steps much smaller than the visible solution’s slow time scale would suggest for accuracy. The handout’s coupled example is

Its distinct real eigenvalues are $(-101\pm\sqrt{9797})/2$, approximately -1.0101 and -99.9899. A diagonalizable modal representation gives slow and fast factors $e^{\lambda_1t}$ and $e^{\lambda_2t}$. The corresponding time constants are about 0.9900 and 0.0100. Requiring both FE modal factors to satisfy $|1+h\lambda_j|\le1$ gives $h\le2/99.9899\approx0.0200$, hence at least 500 steps to reach $t=10$. Scalar eigenvalue tests do not exclude nonnormal transient growth in the original coordinates.

A fast component becoming visually small does not remove its eigenvalue from the equations: roundoff or other perturbations can re-excite it. The reaction example below makes the same restriction observable through accepted steps and work counts.

solve_ivp supplies an ODE interface, $\dot{\mathbf z}=\mathbf f(t,\mathbf z)$. Its BDF option is not the IDAS DAE residual interface. Neither the word “implicit” nor the family name alone establishes the stability of every implementation at every order. BDF1 and BDF2 are A-stable; higher classical BDF orders have smaller stability sectors. Sparse Jacobians, factorization reuse, nonlinear convergence, scaling, and requested accuracy affect practical cost.

Official SciPy solve_ivp documentation, CVODES mathematical description, and IDAS mathematical description describe these interfaces. The work counters below are measured values for the stated experiment, not a universal ranking.

Dense output, events, and checking a run¶

dense_output=True returns an interpolant between accepted mesh points; evaluating sol.sol(t) does not force the solver to step at every requested time. The reaction comparison uses this interpolant on a common checking grid. Interpolation error still depends on the method and tolerances.

An event function returns a scalar whose zeros identify requested occurrences, such as a concentration reaching a threshold. solve_ivp locates sign-changing event roots over accepted steps; several crossings inside one step can be missed. An event may be terminal or restricted by crossing direction. A known control discontinuity should be supplied as an interval boundary and a restart, rather than relying on a smooth-step interpolation across it.

Check success, the diagnostic message, and the reached final time before using a trajectory in optimization. If a terminal event is intended, compare against that event instead of requiring the original final time. Inspect nfev, njev, and nlu together with achieved error, conservation, and the accepted mesh. The official SciPy interface documentation defines these results.

Step changes, control jumps, and method selection¶

An embedded RK pair shares stages between two formulas with different weights. Their difference estimates a local error, which the driver compares with scales such as atol + rtol * abs(y). This does not directly bound the global error at every time. Check important outputs by tightening tolerances, comparing a second method, or using an analytic reference when available.

A one-step RK method requires the current state and its stages. A multistep BDF or Adams method also uses a history. On a nonuniform mesh the coefficients or history representation must change consistently; reusing a constant-step stencil after arbitrarily changing $h$ changes the method.

At a known piecewise-constant control jump $t_s$, integrate exactly to $t_s$, use the prescribed one-sided control on each interval, then restart with the continuous state and the new control. An RHS discontinuity inside a step invalidates the smooth local expansion. Multistep history built across the jump is inappropriate and needs restarting or reducing the method order. A one-step method avoids carrying such history, but still must locate the jump and apply the correct RHS. State-reset events additionally require applying the reset map; ordinary continuous-state control switches do not.

For a small nonstiff ODE with frequent known switches, an explicit RK driver is a reasonable starting point. For a stiff ODE, compare stiff-capable methods with achieved accuracy and Jacobian cost included. For a residual-form DAE, use a suitable DAE solver and consistent initial conditions. For a large PDE discretization, examine sparsity and linear algebra before assuming that fewer RHS calls implies less total work. No single driver is best for all four cases.

RK stages, tableaus, and stability functions¶

A Butcher tableau stores stage dependencies compactly:

The row $a_{kj}$ tells how derivative sample $j$ contributes to stage state $k$; $b_j$ tells how that sample contributes to the endpoint. A strictly lower triangular $\mathbf A$ gives explicit stages in sequence. A lower triangular matrix with nonzero diagonal entries gives implicit stage solves. A full matrix generally couples the stages. Nodes in $[0,1]$ are used by the standard formulas here, but that interval is not a universal definition of RK methods.

FE has $A=[0],b=[1],c=[0]$. BE has $A=[1],b=[1],c=[1]$. Explicit midpoint has

Applying the general RK formula to $z'=\lambda z$ gives stage vector $\hat{\mathbf z}=z_i\mathbf1+\zeta\mathbf A\hat{\mathbf z}$. Where $\mathbf I-\zeta\mathbf A$ is invertible,

This derives the stability function directly from the same coefficients that specify the method. The $h_i$ multiplying every derivative sum is essential for dimensions and consistency. Biegler’s printed equations (9.8) and (9.17), pp. 255 and 259, omit it; the formulas here include it.

Standard Gauss, Radau IIA, and the Lobatto IIIA collocation methods used in the next lecture are implicit RK methods. “All collocation is implicit” would be too broad: a one-node construction at the left endpoint gives FE.

Stiff chemical-reaction example¶

Consider consecutive first-order reactions $A \rightarrow B \rightarrow C$ with time in seconds, $k_1=1000\,\mathrm{s}^{-1}$ and $k_2=1\,\mathrm{s}^{-1}$:

The Jacobian eigenvalues are $-k_1$, $-k_2$ and 0, so the ratio of the two decaying rates is $k_1/k_2 = 1000$: $A$ disappears in about $1/k_1 = 1$ ms while $B$ is still converting a thousand times more slowly. The zero eigenvalue describes conserved total concentration; it is excluded from this ratio. We compare one explicit adaptive method (RK45) with three methods designed for, or able to detect, stiffness (Radau, BDF, LSODA).

Use the same tolerances for every method, and compare solver work through accepted steps and right-hand-side evaluations rather than wall time, which varies substantially across computers and Colab sessions.

All four methods agree on the concentrations to better than 10^-7, but they do not need comparable work. In the archived course run (rerun the preceding cells to obtain counts for your environment):

The mechanism is in the left panel of the figure: once the fast transient is over, the solution is a smooth exponential decay and accuracy alone would permit a large step. Radau, BDF and LSODA take one. RK45 cannot, because its step is capped by the fast eigenvalue $-k_1$ — still -1000 long after the mode it belongs to has died — at the Dormand--Prince stability limit $h \approx 3.3/k_1$. The step it actually takes, $3.25 \times 10^{-3}$, is that ceiling, not an accuracy requirement.

Exact counts depend on software version and tolerances; the comparison is method-specific. RK45 uses about 5.7 times Radau’s right-hand-side evaluations, 19.8 times BDF’s, and 22.8 times LSODA’s in this run. Right-hand-side counts omit Jacobian construction and linear algebra: compare njev, nlu, achieved error, and problem size as well.

For this archived SciPy 1.18.0 run, RK45 makes two initial right-hand-side evaluations and six new evaluations per attempted step. The endpoint derivative is reused at the next attempt (the first-same-as-last, or FSAL, property). Thus $11{,}012=2+6(1{,}835)$ records 1,835 attempts. Of these, 1,582 were accepted and 253 were rejected. The cost is about $11{,}012/1{,}582=6.96$ evaluations per accepted step, rather than exactly six. This accounting uses the stated implementation and automatic initial-step selection; it is not a universal cost formula for every RK driver. See the SciPy 1.18.0 RK implementation.

Discussion: Which method would you try first if this integration were repeated inside an optimization algorithm? What additional information would you check before deciding?

Sources and further reading¶

• Ascher, U. M. & Petzold, L. R. (1998), Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations: pp. 38-39 (normalized local error, zero-stability, convergence), pp. 42-43 and 56-57 (scalar stability and stiffness), pp. 76-80 (quadrature and explicit RK), pp. 102-104 (collocation, stiff accuracy and decay), Exercise 4.13, p. 119 (algebraic stability).

• Biegler, L. T. (2010), Nonlinear Programming: §9.2, pp. 253-260 (integration, RK, multistep methods and DAEs); Theorem 10.1, p. 291 (Gaussian quadrature). Use $h_i$ in each RK derivative sum and differentiate the chosen DAE residual to determine its Jacobian signs.

• Burrage, K. & Butcher, J. C. (1979), “Stability criteria for implicit Runge-Kutta methods,” SIAM Journal on Numerical Analysis 16(1), 46-57, DOI: definitions on pp. 48-49 and Theorem 3.3 on p. 50. The distinct-node condition is needed for the converse equivalence.

• Gear, C. W. & Petzold, L. R. (1984), “ODE methods for the solution of differential/algebraic systems,” SIAM Journal on Numerical Analysis 21(4), 716-728, DOI: Theorem 2.1, p. 719, is the fixed-step BDF convergence result with its regularity and initial-value assumptions.

• Bynum et al. (2021), Pyomo: Optimization Modeling in Python, 3rd ed., §12.1, p. 181: complete ODE-constrained optimization example implemented above.

The derivations distinguish scalar test-equation behavior from nonlinear and DAE claims. The code includes endpoint, residual, conservation, and analytic-error checks; successful integration alone is insufficient evidence of the requested accuracy.