Heat Conduction in a Slab, a Cylinder and a Sphere

Prepared by: Prof. Alexander Dowling (2026)

This notebook is the PDE companion to Lecture 7. We take one transient conduction problem, write it for three geometries at once, discretize space with the method of lines, and solve the resulting differential-algebraic system with Pyomo.dae.

The physical problem and the idea of indexing the three geometries with a single shape parameter are adapted from J. C. Kantor, ND Pyomo Cookbook, notebook 5.03 “Heat Conduction in Various Geometries”. The derivation, the code and the numerical experiments below are our own.

Learning Objectives¶

• Nondimensionalize a transient conduction PDE and recognize the Fourier number as its only time scale

• Write one model that covers a slab, a cylinder and a sphere through a shape index $s$

• Apply the method of lines: discretize space first, then hand a DAE in time to Pyomo.dae

• Compare finite difference, finite volume and finite element spatial discretizations, and say where each belongs

• Handle the $r = 0$ coordinate singularity correctly

• Run the IDAES Diagnostics Toolbox on a Pyomo.dae model and read its report -- including the difference between a warning and a caution

• Recognize a dangling variable: a derivative that appears in no discretization equation, which IPOPT will happily declare optimal

Use the boundary-condition guide for Dirichlet, Neumann and Robin equations, the center limit for regularization, and the boundary-control extension for stiffness and active path constraints. The executable FD and FV builders below are the correct models; the alternative-discretization tables retain historical failure evidence, without adding runnable broken variants.

The Physical Problem¶

A solid body, initially at uniform temperature $T_0$, is dropped at $t=0$ into a bath held at $T_\infty$. The surface reaches the bath temperature immediately, and we want the temperature everywhere inside as a function of time.

With constant positive thermal conductivity $k$, density $\rho$ and heat capacity $C_p$, no source term, and no axial or angular temperature variation, the energy balance is

Divide by $\rho C_p$ and define the thermal diffusivity $\alpha = k / (\rho C_p)$, which has units of length$^2$/time:

One Equation, Three Geometries¶

The three classical bodies -- an infinite slab of half-thickness $R$, an infinite cylinder of radius $R$, and a sphere of radius $R$ -- are all one-dimensional: temperature depends on a single space coordinate $r$ measured from the center. In each case the Laplacian of a radially symmetric field is

where the shape index $s$ is

For the slab, $s=0$ and the extra term vanishes -- Cartesian coordinates. For $s = 1$ and $s = 2$ the geometry contributes the $\tfrac{s}{r}\partial T/\partial r$ term, which is singular at the center. That singularity is the interesting part of this notebook, so hold onto it.

Nondimensionalization¶

Assume $R>0$ and $T_\infty\ne T_0$. If the two temperatures agree, this normalization is undefined and the constant-temperature solution should be treated directly. Three scales are available: the temperature difference the body must cross, the body’s own size, and the time the heat needs to cross it. Define

Substituting, with $T = T_0 + (T_\infty - T_0) T'$, $r = R r'$ and $t = R^2 t' / \alpha$:

Every term carries the same factor $\alpha (T_\infty - T_0)/R^2$, which therefore cancels:

with

Two things are worth noticing.

1. No parameters survive. The dimensionless problem depends only on $s$. Solve it once per geometry and rescale it for other constant positive material properties, sizes and distinct temperature pairs under these assumptions. Variable properties, sources or a convective boundary introduce additional parameters.

2. $t'$ is the Fourier number, $\mathrm{Fo} = \alpha t / R^2$: the ratio of the time elapsed to the time heat needs to diffuse across the body. $\mathrm{Fo}$ of order one means diffusion has acted over the body scale; a specified closeness to equilibrium still depends on the geometry and the tolerance.

The last condition is the symmetry (or no-flux) condition. It is not an extra physical assumption -- it follows from the body being symmetric about $r'=0$, so no heat can cross the center. It is, however, the boundary condition that will cause all the trouble.

From here on we drop the primes and write $T$, $r$, $t$ for the dimensionless quantities.

Dimensionally, $[\alpha]=\mathrm{length}^2/\mathrm{time}$, so $\alpha t/R^2$ is dimensionless, whereas $t/R$ and $\alpha t$ are not. A diffusive penetration length scales as $\sqrt{\alpha t}$; the Fourier number is its squared ratio to $R$. This is a scale estimate, not a sharp front location.

The Method of Lines¶

We have one equation in two independent variables. Pyomo.dae is comfortable with either, but the standard strategy for a parabolic PDE is to break the symmetry between them:

Method of lines. Discretize the spatial derivatives on a fixed mesh $0 = r_0 < r_1 < \dots < r_N = 1$. What is left is a system of ODEs (here, a DAE, because the boundary rows are algebraic) in time alone, one equation per mesh point -- the “lines” of the name.

Writing $T_i(t) \approx T(t, r_i)$, the semi-discrete system is

with $\mathcal{L}_i$ the discrete approximation to $\partial^2 T/\partial r^2 + (s/r)\,\partial T/\partial r$ at node $i$, plus algebraic rows at $i = 0$ (symmetry) and $i = N$ (the bath temperature).

Why this order?

• The time direction is an initial-value problem, and Lecture 7’s machinery -- orthogonal collocation on finite elements -- is built exactly for that. We can hand the time direction to dae.collocation with right Radau. Its ODE stability properties are useful for diffusion, but temporal accuracy, algebraic consistency, nonlinear convergence and spatial error still need checking.

• The space direction is a boundary-value problem in disguise. All of its structure -- the boundary conditions, the coordinate singularity, the conservation properties -- is decided by how we build $\mathcal{L}_i$. Doing it deliberately, rather than letting a transformation do it for us, is the whole game.

In Pyomo.dae both directions are ContinuousSets and both get a TransformationFactory, so the method of lines is a way of thinking about the model rather than a different API. But it tells you which transformation to worry about.

The reduction does not guarantee a regular index-1 DAE. Algebraic boundary equations or active constraints can raise the index. The smooth-ODE endpoint orders apply only under the relevant regularity and consistency assumptions. Ascher & Petzold, Example 1.7, pp. 12–13, explains why a spatially fixed method of lines is especially useful when time and space play different roles; rapidly moving sharp fronts may require another mesh strategy.

Three Ways to Discretize Space¶

For this linear conduction model, all three families give a sparse spatial operator on discrete unknowns (nodal values or cell values). They differ in what they ask the approximation to reproduce exactly, and that choice is what makes one appropriate and another awkward for a given problem.

Finite difference (FD): match the Taylor series¶

Pick a stencil and choose coefficients so the truncation error is as high order as possible. On a uniform mesh of spacing $h$,

• Cheap and transparent. Every row is three numbers. This is what dae.finite_difference builds.

• Appropriate when the geometry is simple, the mesh is structured, and the coefficients ($k$, $\rho C_p$) are smooth.

• Awkward when you need a guarantee. FD approximates the differential equation pointwise; pointwise consistency alone does not establish a discrete conservation law. Conservative FD formulations do exist, and simple FD/FV stencils can coincide. At a coordinate singularity use a regularized equation or an integral balance, rather than evaluating $1/r$ at zero.

Finite volume (FV): match the integral balance¶

Partition the domain into control volumes (here, shells) and integrate the conservation law over each one. Using the self-adjoint form $\partial_t T = r^{-s}\,\partial_r ( r^{s} \partial_r T )$ and integrating $r^{s}\,dr$ over the cell $[r_{i-1/2},\,r_{i+1/2}]$,

Approximating each face gradient by a two-point difference gives a scheme in which the flux leaving one cell is identically the flux entering the next, so the sum over cells telescopes: the semi-discrete balance conserves energy algebraically on any mesh. Numerical time integration and solver residuals determine how closely a computed trajectory satisfies that balance.

• Appropriate when conservation matters (reactors, flow, shocks), the mesh is unstructured, or the coefficients jump at an interface -- FV puts the interface on a face, where the flux is single valued.

• Note what it does at $r = 0$: the inner face of the first cell has area $r_{1/2}^{\,s} \to 0^{s} = 0$ for $s > 0$. The symmetry condition is structural: no special case, no boundary condition to write. We will use this below.

Finite element (FE): match a weighted residual¶

Expand $T(t,r) \approx \sum_j T_j(t)\,\phi_j(r)$ in basis functions and require the residual to be orthogonal to a set of test functions. Galerkin FE integrates the weak form by parts, which lowers the smoothness required of the trial functions. For self-adjoint diffusion with appropriate boundary conditions and quadrature, the resulting stiffness matrix is symmetric and an energy argument establishes stability. These properties are not universal for every FE model. Orthogonal collocation instead sets the residual to zero at selected points inside each element; it is not the same discretization as Galerkin FE.

• Appropriate when the geometry is complicated, the solution is smooth enough to reward high order, or you want error estimates and adaptivity on a rigorous footing.

• In Pyomo.dae, dae.collocation supplies the collocation construction, not a general Galerkin FE assembler. We use collocation in time below. It can also be applied in space -- and, as we will see, it inherits exactly the same $r=0$ trap as the finite difference scheme.

A shared residual viewpoint¶

For a trial field $T_h$, let $\mathcal R(T_h)$ be the differential-equation residual. A weighted-residual method asks for $\int W_j\mathcal R(T_h)\,dr=0$ for each test weight. Galerkin uses basis/test functions and typically integrates a diffusion term by parts. FV uses one subdomain indicator at a time and evaluates face fluxes, giving the control-volume balance. Collocation samples the strong residual at points, formally corresponding to point-evaluation (Dirac) weights. FD instead derives local derivative approximations from Taylor expansions. These viewpoints can yield the same stencil in a special case, but they impose different conditions in general. See Patankar, Section 3.2, pp. 27–31.

The Coordinate Singularity at $r = 0$¶

For $s > 0$ the term $\dfrac{s}{r}\dfrac{\partial T}{\partial r}$ is undefined at $r = 0$. The solution is perfectly well behaved there -- by symmetry $\partial T/\partial r \to 0$ at least as fast as $r$, so L’Hopital’s rule gives $\lim_{r \to 0} \frac{s}{r}\frac{\partial T}{\partial r} = s\,\frac{\partial^2 T}{\partial r^2}$ -- but the equation as written cannot be evaluated at $r=0$. Two consequences for the discrete model:

1. Do not evaluate the unregularized $1/r$ expression at zero. For a smooth symmetric solution, $T_r(t,0)=0$ and $T_r(t,r)=rT_{rr}(t,0)+o(r)$. Hence the finite center equation is

   $$T_t(t,0)=(1+s)T_{rr}(t,0).$$

   (11)

   One valid discretization uses that regularized evolution equation with a symmetry-aware stencil. The builder below instead skips the center PDE and supplies a second-order Neumann row tying the center temperature to nearby temperatures. These are different closures; do not add both indiscriminately. For the skipped-row formulation, the derivative must have a defining equation as well as the prescribed value.

2. Skip the corner $(t, r) = (0, 1)$. The initial condition asks for $T = 0$ there and the boundary condition asks for $T = 1$. Both cannot hold. We keep the boundary condition (the bath is imposed from $t=0^+$ onward) and skip the initial condition at $r=1$. Physically this is the discontinuity that makes the surface flux infinite at $t = 0^+$; numerically it is simply one row we must not write twice.

The second point is routine bookkeeping. The first one hides a trap.

The builder uses a central interior stencil and a one-sided second-order center row. Read the defining equation after the transformations as well as the physical boundary condition; both are needed for this choice of closure.

A Pyomo.dae Model¶

Here is the model, written once for all three geometries. Read the three constraint rules first -- they are a direct transcription of the boxed PDE and its three side conditions. Then read the block after the transformations, which supplies the one equation the spatial discretization does not write for us. The next section explains why it is needed and what happens without it.

The Trap: A Derivative That Appears in No Equation¶

Everything in build_fd up to the transformations is a transcription of the mathematics. The block after them is not, and it is the part that is easy to leave out. Here is the natural way to write the symmetry condition:

@m.Constraint(m.t)
def bc_center(m, t):
    return m.dTdr[t, 0] == 0        # DO NOT DO THIS AND NOTHING ELSE

It reads exactly like $\partial T/\partial r|_{r=0} = 0$. It is also wrong, for a reason that has nothing to do with the mathematics and everything to do with what scheme="CENTRAL" actually generates.

A central difference needs a neighbour on both sides. At an interior node it has one:

At $i = 0$ there is no $T_{-1}$, so Pyomo.dae writes no equation at all for $\texttt{dTdr}[t, 0]$. It is not set to a one-sided difference; it is not set to anything. It is a variable that appears in no discretization equation -- a dangling variable. The constraint above would then not impose symmetry on the temperature field; it would assign a value to a number that touches nothing else in the model, and the center temperature would be determined by nothing at all.

We can see the hole directly, without building a model at all. Discretize a bare DerivativeVar and ask which mesh points the generated defining equations cover.

Read the CENTRAL row: neither endpoint gets a defining equation for the first derivative, and neither gets one for the second. At $r=1$ that is harmless -- the PDE row is skipped there and the bath temperature is imposed directly. At $r=0$ it is fatal, because the PDE row is skipped there too.

The other two rows are the reason “just use a one-sided scheme everywhere” is not a fix: FORWARD and BACKWARD each close one hole and open another somewhere else. They move the problem rather than solve it. See the first table below.

So build_fd writes both halves of the symmetry condition itself: it fixes $\texttt{dTdr}[t,0]$ to zero, which is the physics, and it adds bc_center_def, which is the equation the scheme never wrote. Together with $\texttt{dTdr}[t,0] = 0$ the three-point stencil collapses to $3T_0 = 4T_1 - T_2$, the second-order statement that the profile is flat at the center.

The Reference Solution¶

Each of the three geometries has a classical separation-of-variables solution, so we can measure error rather than eyeball it. For the dimensionless problem above:

Checking the Model Before You Trust It¶

This is a simulation: the temperature field is completely determined by the data, so a correctly posed model must have zero degrees of freedom. That is the cheapest check available, and it happens before the solver is ever called.

The IDAES Diagnostics Toolbox¶

Counting degrees of freedom tells you how many variables are unaccounted for. It does not tell you which. On a model with tens of thousands of rows that difference is the whole game, so run the structural analysis rather than the arithmetic.

The IDAES Diagnostics Toolbox does exactly this. report_structural_issues() needs no solution and no initialization -- it works on the incidence pattern of the model, so you can run it the moment the model is built. Make it the first thing you do to any Pyomo.dae model, before you ever look at a number.

Zero warnings. No structural singularity, no under-constrained set, no over-constrained set, degrees of freedom consistent with the constraint count. That is what a well-posed Pyomo.dae model looks like.

The two cautions are worth understanding, because they are not defects and you will see them on every model of this shape.

• 13 variables fixed to 0 -- these are $\texttt{dTdr}[t,0]$ at the 13 time points, fixed by build_fd as statement (1) of the symmetry condition. The toolbox flags variables fixed to zero because a stray fix(0) is a common way to accidentally delete physics. Here it is the physics.

• 49 unused variables -- variables that appear in no activated constraint. All 49 are inert: 9 are $\texttt{dTdt}$ at $t = 0$, which Radau collocation never references; 26 are $\texttt{d2Tdr2}$ on the two boundary rows $r = 0$ and $r = 1$, where the PDE is skipped; and 14 are $\texttt{dTdr}$ at the surface and at $(t,r) = (0,0)$, likewise unreferenced. Pyomo.dae creates a full grid of derivative variables and only some of them are used. They influence no result.

A caution is an invitation to look, not a finding. Knowing which cautions your model should produce is what lets you notice the day it produces a different one.

Alternative Discretizations at $r = 0$, and What They Do¶

The table below documents the discretizations we did not ship, so that the failure modes are on record without a broken model being a copy-paste away.

All numbers were measured on the slab ($s = 0$) unless stated, at $\mathrm{Fo} = t_f = 0.5$, with nfe_t=30, ncp=3 Radau collocation in time, and IPOPT 3.13.2 with ma27. “Center error” is $|T(t_f, 0) - T_\text{exact}(t_f, 0)|$ with $T_\text{exact}(0.5, 0) = 0.6292$; the five values run over $h = 1/10,\, 1/20,\, 1/40,\, 1/80,\, 1/160$ and “order” is $\log_2$ of the ratio of the last two.

The Answer Does Not Depend on the Initial Guess¶

A model with unconstrained degrees of freedom and no objective has no unique answer, and the quickest way to see that is to change the initial guess and re-solve. Run the same check on the correct model: if the answer moves, something is not pinned.

Identical to every digit printed across a range of guesses that includes a physically impossible one, in every geometry and on every mesh. The remaining discrepancy against exact is truncation error, and unlike a structural defect it goes away under refinement -- which is what the next section measures.

Mesh refinement checks a different error¶

Consistency alone does not prove convergence: stability, well-posed boundary equations, a converged temporal discretization and sufficiently small nonlinear-solver error are also needed. Refine space while holding the temporal error small, then refine time on a fixed spatial mesh. Compare against an independent reference, not only the previous mesh.

A flat error curve can indicate a structural defect, but it can also be a temporal- or solver-error floor. Inspect the equations and vary those tolerances before deciding which error dominates.

Read the order column: it is $\log_2$ of the ratio of the last two errors, so 2.0 means second order in $h$, 1.0 means first order, and $\approx 0$ would mean the error is not going anywhere. All three geometries land at 1.9--2.0, which is what the interior stencil and the center row were each built to deliver -- and it is only true because they deliver it at the same order. The first-order center row in the table above produces the same picture with slope 1, and the dangling model produces a flat line.

The curves approach second-order behavior on the finer meshes; the reported order uses only the last pair. Coarser pairs, particularly for the slab, can be pre-asymptotic. Read the actual pairwise ratios rather than interpreting every segment as a line of slope 2. An observed convergence rate is evidence for this experiment, not a general proof of the scheme.

Reproducibility and archived figure data¶

The canonical convergence and profile figures are generated from committed numerical data using the plotting cells below. The handout can reuse either figure without introducing another implementation of the PDE. figures/results/heat-convergence.json stores the raw mesh- refinement errors for all three geometries, the initial-guess experiment, and the expected structural-diagnostics summary. figures/results/heat-slab-profile.json stores every profile and analytical-series comparison for the slab, cylinder, and sphere; the dedicated profile figure uses the slab subset.

A maintainer re-renders either handout figure without Pyomo or IPOPT with:

python figures/render_from_notebook.py heat-convergence heat-slab-profile

The archived runs used Pyomo 6.10.1, IPOPT 3.13.2 (IDAES build with ma27), nfe_t=30, and three Radau collocation points for the refinement experiment; the profile data use nfe_r=80, nfe_t=50, and three Radau points. Exact iteration counts and timings are intentionally not used as scientific evidence: they depend on the operating system, solver build, linear solver, and CPU. The verification targets are: the degree-of-freedom count, expected diagnostics cautions, initial- guess invariance, and second-order convergence are the checks to reproduce.

The Correct Model, Geometry by Geometry¶

The builder is the same in all three cases. What changes is why each piece is correct.

Slab, $s = 0$¶

The PDE is $T_t=T_{rr}$ without a coordinate singularity, but the half-slab still requires symmetry at its midplane. A center evolution row is possible only with a stencil that encodes that symmetry, for example by an even ghost value. Merely dropping the Neumann information and writing the PDE does not close the problem. The shipped model instead uses one center Neumann row with its defining stencil, and skips the center PDE, for all three geometries.

Cylinder, $s = 1$¶

The unregularized $\frac{1}{r}\partial_r T$ cannot be evaluated at zero. Skip it as this builder does, or replace it by the regularized center equation derived above. The symmetry condition is doing real work -- it is the only equation that determines the axis temperature, and without a defining equation for $\texttt{dTdr}[t,0]$ the axis is unconstrained. Every interior row $i = 1, \dots, N-1$ has $r_i = i h > 0$, so the $1/r$ coefficient is finite everywhere it is used; it is large ($1/h$) at the first interior node, which is why the discrete operator is stiff near the axis and why an $\mathcal{O}(h^2)$ symmetry stencil is worth having.

Sphere, $s = 2$¶

Identical structure, with $2/r$. The $1/r$ weighting is stronger, so the near-center rows dominate the Jacobian even more -- which is also why the dangling model of the first table is least visibly wrong here. Nothing about the correct model changes.

In all three cases the count works out the same way, per time point: $N$ unknowns $T_0, \dots, T_{N-1}$ (the surface value $T_N$ is imposed), $N-1$ PDE rows at the interior nodes, and one center row. $N$ equations, $N$ unknowns, zero degrees of freedom -- which is what the DOF table above confirms.

The dedicated slab profile below shares the archive used for the three-panel comparison. The slab is especially useful for exposing the documented missing-center-row defect; cylinder and sphere results remain available above. Both the full profile and convergence experiments belong in this reference guide even when lecture space permits only one figure.

At a fixed Fourier number the sphere is furthest along and the slab is furthest behind. That is the surface-to-volume ratio: per unit volume, a sphere presents $3/R$ of surface, a cylinder $2/R$, a slab $1/R$. The shape index $s$ is doing exactly that bookkeeping, and it is the reason $\mathrm{Fo}$ alone is not enough to say “how hot is it inside” without also saying which shape you meant.

The Finite Volume Alternative¶

Everything above was finite difference. It is worth building the finite volume version of the same problem, because it makes a structural point that no amount of care with stencils can:

In the finite volume formulation the symmetry condition at $r = 0$ is not a boundary condition at all. It is a consequence of the geometry, and there is nothing to forget.

Place $N$ shells with faces at $r_{i+1/2} = (i+1)h$, $h = 1/N$, and unknowns at the cell centers. The balance on shell $i$ is the discrete form of the integral law from earlier,

For the innermost shell, $r_{-1/2} = 0$ and therefore $A_{-1/2} = 0^{s} = 0$ for $s > 0$: the inner face has no area, so no heat can cross it. For the slab, $s = 0$ makes $A_{-1/2} = 1$, and we set the inner flux to zero by symmetry -- one line, and the only geometry-specific line in the model.

At the outer face of the last shell the temperature is known, and the half-cell distance from the center of that shell to the surface is $h/2$, so the flux is $A_{N-1/2}\,(1 - T_{N-1}) / (h/2)$.

Note also that there is no DerivativeVar in space at all. The method of lines is explicit here: the only ContinuousSet is time.

Zero degrees of freedom in every geometry, second-order convergence, and no symmetry constraint anywhere in the model. The Diagnostics Toolbox says the same thing structurally.

The FV diagnostics again report zero structural warnings. For the cylinder and sphere, bounded center gradients multiply a zero-area inner face, making its net heat flux zero. For the slab, the face has nonzero area and flux_in = 0 explicitly imposes symmetry. Finite volume does not prevent an incorrect boundary flux or an incorrect geometry model; it makes the balance and the boundary-flux choice visible.

Cell-centered FV has no center temperature variable requiring a separate Neumann derivative row. Its unknowns are cell values rather than boundary values; surface reporting and boundary fluxes therefore need the half-cell distance. The following table retains historical measured failures of alternative formulas. Those variants are not runnable builders in this guide.

Alternative Finite Volume Discretizations, and What They Do¶

Like the first failure table, these are historical measurements under the stated configuration, not fresh executions in this revision; none ships as a runnable option above. Conditions are $\mathrm{Fo} = t_f = 0.5$, nfe_t=30, ncp=3 Radau, IPOPT 3.13.2 with ma27, and $N = 10,\, 20,\, 40,\, 80$ cells. The error is taken at the first cell center, $r = h/2$, against the series solution evaluated there; “order” is $\log_2$ of the ratio of the last two.

What each boundary condition supplies¶

Dirichlet prescribes the value, such as $T(t,1)=1$. At that boundary, the value is already supplied; the interior PDE stencil need not create a second evolution equation for it. For a boundary temperature control $u(t)$, write $T(t,1)=u(t)$ and include the control’s bounds, cost and parameterization.

Neumann prescribes a normal derivative or flux, for example $T_r(t,0)=0$. It does not directly set the temperature. If the PDE row is omitted, a defining stencil must tie that derivative to the temperature unknowns. The shipped center row is

so $T_0=(4T_1-T_2)/3$. This is why fixing a dangling DerivativeVar to zero alone fails. Skip the same row at the initial time if the initial conditions already imply it; adding dependent equations can destroy a square solve.

Robin couples temperature and outward heat flux. Write the physical convective boundary as

After the stated scaling,

Here $h_c$ is a heat-transfer coefficient, distinct from a time-element length. With a cold solid and hot bath, $T_r(t,1)>0$, so the outward heat flux is negative: heat enters the solid. A central stencil also leaves the derivative at the outer boundary undefined. Supply a one-sided defining row or a boundary-face flux derived from the balance; the Robin equation alone does not create that row. The large-Biot limit approaches the fixed-surface-temperature model under the usual regularity assumptions.

Row rule. Omit a PDE or derivative row only when another equation supplies the missing information. Count unknowns and independent equations after transformation, inspect the derivative coverage, then check the residual and physical behavior. A zero count alone cannot detect dependent equations. The IDAES diagnostics workflow separates structural checks from checks requiring a numerical solution.

Initial-boundary corner. The interior initial data are $T(0,r)=0$ for $r<1$, while the boundary jumps to 1 for positive time. They are incompatible at $(0,1)$. The code retains the boundary value there and skips the conflicting initial row. It does not assert a smooth solution at that corner; early boundary layers need a finer or graded spatial mesh and temporal resolution.

Boundary controls, spatial stiffness and active path constraints¶

For Biegler Example 10.6 (pp. 318–319), let $z_t=z_{xx}$ on $0<x<\pi$, $0<t\le5$, with $z(x,0)=0$ in the interior and boundary controls $z(0,t)=u_0(t)$, $z(\pi,t)=u_\pi(t)$. The cost penalizes field temperature and controls:

A path lower bound is $z(x,t)\ge\sin(x)\sin(\pi t/5)-0.7$. This is a distinct optimal-control extension, not a property of the fixed-bath simulation.

With $x_j=j\delta$, $\delta=\pi/n$, the interior ODE is

The matrix is tridiagonal; only the first and last rows directly see the boundary controls. Its homogeneous eigenvalues are

The fastest decay grows like $\delta^{-2}$ while the slowest approaches a finite rate: spatial refinement increases stiffness. For forward Euler on this uniform diffusion system, $\Delta t\le\delta^2/2$ is a sufficient absolute-stability restriction. Right-Radau time integration removes that explicit stability restriction for this linear system, but accuracy and algebraic consistency still matter.

To match Biegler’s scaled time, set $\tau=t/\delta^2$; the ODE loses its $1/\delta^2$ factor and the horizon becomes $5/\delta^2$. Trapezoidal spatial quadrature gives the scaled cost

The bound becomes $z_j(\tau)\ge\sin(j\delta)\sin(\pi\delta^2\tau/5)-0.7$, with endpoint lower bounds $u_0,u_\pi\ge-0.7$. Both the time differential and the explicit time argument change; scaling only the ODE is incomplete.

If an interior path constraint is active over an interval, treating it as an equality can create a high-index constrained system. Repeated differentiation propagates information from that spatial node toward the boundary controls; for the symmetric center constraint in this example the number of differentiations grows with spatial resolution (Biegler p. 319 gives $n/2$). This is why a general promise of smooth-ODE order for every discretized PDE is too strong. The fixed-bath simulations here contain no such active path constraint.

For the method of lines and eigenvalues, see LeVeque (2007), Sections 9.2–9.4, pp. 184–186; for boundary control and index growth, Biegler (2010), Example 10.6, pp. 318–319.

Check time integration against the semi-discrete matrix exponential¶

The analytic series checks the continuous PDE. A second reference isolates the temporal error: eliminate the center temperature using $T_0=(4T_1-T_2)/3$, form the interior FD matrix directly, and solve this linear ODE with a matrix exponential. The equilibrium is the all-ones vector, so $\mathbf T(t)=\mathbf1-\exp(t\mathbf L)\mathbf1$. This calculation uses no Pyomo derivative transformation and no collocation coefficients.

Agreement between the matrix exponential and the transcribed trajectory checks the time solution of this spatial model. Their common difference from the series is a spatial error. Neither reference should be used to relabel a solver’s successful termination as a proof of PDE correctness.

Spatial-method references¶

LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations (2007), Chapter 2, pp. 20–21, states convergence in terms of consistency and stability; Chapter 9, pp. 184–186, develops the diffusion method of lines and its stiffness. Patankar, Numerical Heat Transfer and Fluid Flow (1980), Section 3.2-4, pp. 30–31, develops the control-volume balance and its conservation property. The weighted-residual viewpoint explains Galerkin FE, FV subdomain integration, and collocation as distinct residual choices; it does not make their stencils interchangeable.

The full dynamic-transcription derivation is in Dynamic Optimization with Collocation and Pyomo.DAE. Exact collocation coefficients and separate stability tests are in Collocation Points and Runge-Kutta Coefficients.

Take-Home Messages¶

1. Nondimensionalize first. Under the constant-property, symmetric, fixed-surface-temperature assumptions, the three geometries reduce to one equation in one parameter $s$, with the Fourier number as the clock.

2. The method of lines is a decision, not a syntax. Discretizing space first is what lets you hand the time direction to Lecture 7’s collocation machinery and concentrate on the spatial operator, where the physics lives.

3. Finite difference, finite volume and finite element differ in what they reproduce exactly -- the Taylor series, the integral balance, or a weighted residual. Reach for FV when conservation or an interface matters, FE when the geometry or the accuracy demands it, FD when neither does.

4. 🔴 A Neumann boundary condition written on a DerivativeVar may constrain nothing. CENTRAL finite differences and Radau collocation both leave the derivative at the left endpoint of the domain with no defining equation. Writing m.dTdr[t, 0] == 0 and stopping there assigns a value to a dangling variable, and the center temperature is unconstrained. Write the defining equation too.

5. optimal is not correct. Every failure in the two tables above terminated cleanly, several of them with plausible-looking numbers. Three things catch this class of defect: a degree-of-freedom count -- a regular determined simulation should have zero -- report_structural_issues(), which says which variables are loose rather than how many, and mesh refinement, which turns a structural defect into a visible plateau. Do all three.

6. Run the diagnostics before you run the solver. report_structural_issues() needs no solution and no initialization; it costs seconds and it reads the incidence pattern you cannot read by eye. Learn what cautions your model should emit, so that a new one means something.

7. Test on the geometry where the bug is loud. The finite difference defect was obvious in the slab and nearly invisible in the sphere; the finite volume defect is the other way round. A model verified only on the forgiving case is not verified.

Try It Yourself¶

1. Predict, before running it, which endpoint loses its defining equation under scheme="BACKWARD" and what that breaks in build_fd. Check your answer with stencil_coverage("BACKWARD") -- and then say what one-sided row you would owe the model to repair it.

2. Replace dae.finite_difference in space with dae.collocation in a copy of build_fd and confirm, with the toolbox rather than by eye, that the derivative at $r=0$ still needs a defining equation of your own.

3. Add a Robin (convective) boundary condition at the surface, $-\partial_r T|_{r=1} = \mathrm{Bi}\,(T(t,1) - 1)$, where $\mathrm{Bi} = h_cR/k$ is the Biot number. Which endpoint does CENTRAL leave undefined now, and what does that imply about which one-sided stencil you owe the model? Recover the results above as $\mathrm{Bi} \to \infty$.

4. Grade the mesh toward $r=1$ so the short-time boundary layer is resolved, and redo the $\mathrm{Fo}=0.01$ comparison. ContinuousSet(initialize=[...]) accepts an explicit list of points.

5. Reproduce one row of either table for yourself in a scratch copy of the builder -- then delete it. Confirming a documented failure mode is worth an hour; leaving a runnable broken model in a notebook someone else will copy is not.