Stiff CSTR Integration: Stability Step Versus Accuracy Target
A lab-scale CSTR is modelled by dy/dt = -80(y - 1.02) + 0.4 sin(2t), y(0) = 0, where y is dimensionless concentration and t is in minutes. A student integrates with explicit Euler using h = 0.02 min and obtains a smooth response. A colleague uses h = 0.05 min and sees a growing sawtooth oscillation before the computation overflows. The supervisor wants a production run to t = 20 min with an error below 1e-2 in y and suggests either halving h or switching to implicit Euler. The nonlinear term is Lipschitz with constant L = 80. The cost per step for implicit Euler is about five times that of explicit Euler, and the existing code has no Newton solver. The team must choose a method and step size before the next shift.
1. Using the linear test equation, determine the largest stable step for explicit Euler on the fast mode. Explain why h = 0.02 is smooth while h = 0.05 diverges.
Read the worked answer
The fast mode corresponds to the eigenvalue λ = -80. For explicit Euler on y' = λy, the amplification factor is R = 1 + hλ, and stability requires |1 - 80h| ≤ 1. This gives -1 ≤ 1 - 80h ≤ 1, hence 0 < h ≤ 2/80 = 0.025 min. At h = 0.02, R = 1 - 1.6 = -0.6, so the fast component alternates sign but decays by a factor 0.6 per step; the slow forcing dominates and the plotted response looks smooth. At h = 0.05, R = 1 - 4 = -3, so the fast component grows by a factor 3 each step with alternating sign. That is the observed sawtooth, and it eventually overflows. The nonlinearity is Lipschitz with L = 80, so the same estimate is a reliable local stability guide. Note that stability is not the same as accuracy: even at h = 0.02 the first-order truncation error must still be checked against the 1e-2 target.
2. Estimate the number of steps and relative computational cost for (a) explicit Euler with h = 0.02, (b) implicit Euler with h = 0.05, and (c) implicit Euler with h = 0.1. Which would you recommend if the error target is strict?
Read the worked answer
For t = 20 min: (a) explicit with h = 0.02 requires 1000 steps; at unit cost per step the cost is 1000. (b) implicit with h = 0.05 requires 400 steps; at five times the cost per step the equivalent cost is 2000. (c) implicit with h = 0.1 requires 200 steps; equivalent cost is 1000. The implicit method is A-stable, so all three implicit steps are stable in principle. However, implicit Euler is only first order, and at h = 0.1 the local error is O(h^2) with a larger constant, while the accumulated error over 20 min may exceed 1e-2. There is also a solver issue: for h = 0.1, hλ = -8, so fixed-point iteration diverges and a Newton solver would be needed. With the existing code and a strict error target, explicit Euler at h = 0.02 is the lowest-risk choice; if a Newton solver is added, implicit Euler at h = 0.02 or a higher-order implicit method would be preferable.
3. The supervisor says 'stiffness just means use a small step'. Critique this statement using the eigenvalue and the cost numbers.
Read the worked answer
The statement confuses stability with accuracy. The eigenvalue -80 forces explicit Euler to use h ≤ 0.025 min even though the slow forced response has a period of about π min (from sin 2t). The fast mode decays in a time scale of 1/80 = 0.0125 min, but it still dictates the explicit step for the whole 20 min run. That is the practical meaning of stiffness: a stable but fast mode controls the step size, not the accuracy of the slow solution. Implicit Euler removes the stability restriction, which is why h = 0.05 or 0.1 can be considered at all. But the cost numbers show that a five-fold cost per step can erase the benefit: 400 implicit steps at five units each cost 2000, twice the 1000 units for 1000 explicit steps. The correct comparison is stability-limited explicit cost versus accuracy-limited implicit cost, and the latter needs a Newton solver and error control before it can be recommended.