Investigating runge-kutta methods for stiff differential systems with applications to chao phraya urban heat island effects

Introduction

What happens when a well-known numerical recipe meets a problem it wasn't designed for? Explicit Runge-Kutta integrators sit at the heart of most introductory courses on differential equations, yet they falter once the eigenvalue spread of the Jacobian becomes large enough to force impractically small time steps. Stiff systems of this kind arise naturally in reaction kinetics, circuit simulation, and, as this research highlights, in the coupled heat-transfer equations governing urban microclimates. The urban heat island (UHI) effect, the measurable temperature surplus a city maintains over its rural surroundings, has been documented in Bangkok's Chao Phraya corridor since the late 1990s. Surface energy budgets in that region involve fast radiative exchanges layered on top of slow conductive processes in concrete and asphalt, producing a stiffness ratio that can exceed 10⁴ during clear-sky nights [1] . Modelling this behaviour with ordinary differential equations (ODEs) obtained after spatial semidiscretisation demands solvers that remain stable without being forced into tiny time increments. Several implicit Runge-Kutta families address stiffness directly. Diagonally implicit methods (SDIRK) reduce each stage to a single nonlinear solve [2] , while fully implicit Radau and Gauss-Legendre schemes offer higher stage order at greater algebraic cost [3] . Composite methods such as TR-BDF2 blend a trapezoidal half-step with a backwarddifference corrector, achieving L-stability with moderate implementation effort [4] . Despite this variety, head-to-head comparisons on real geophysical stiff systems remain scarce. Earlier work by Tanaka and Kittipong [5] applied a single SDIRK scheme to a simplified one-dimensional UHI model but did not vary solver families. Montoya [6] benchmarked implicit methods on chemical kinetics stiffness but used synthetic eigenvalue spectra rather than physically derived systems. And while Pham et al. [7] surveyed Bangkok's heat island with a mesoscale atmospheric code, they employed a built-in solver without isolating the ODE integration component. This research fills that gap. We formulate a semi-discrete energy-balance ODE system from station data, apply five Runge-Kutta variants spanning the explicit-to-fully-implicit spectrum, and evaluate them on accuracy, stability, and computational cost. The objectives are: (a) to quantify which solver family best matches the stiffness profile of Chao Phraya UHI equations, and (b) to provide practical guidance for selecting integrators in similar urban-climate ODE models.

Theoretical background

The semi-discrete system takes the form dy/dt = f(t, y) where y ∈ ℝⁿ collects nodal temperatures and f encodes diffusive fluxes plus radiative and anthropogenic source terms. After spatial discretisation on a uniform mesh of spacing Δx, the Jacobian ∂f/∂y becomes a banded matrix whose eigenvalues λ satisfy Re(λ) < 0 for all physically meaningful parameter sets [8] . The stiffness ratio S = |λ_max|/|λ_min| exceeded 1.3 × 10⁴ in our Bangkok configuration during nocturnal clear-sky episodes. An s-stage Runge-Kutta method is characterised by its Butcher tableau (A, b, c). When A is lower-triangular the method is diagonally implicit (DIRK); when A is full the method is fully implicit. A-stability, meaning the numerical solution decays whenever Re(λ) < 0, requires that the stability function R(z) satisfy |R(z)| ≤ 1 for all z with Re(z) ≤ 0 [9] . Lstability adds the condition |R(z)| → 0 as Re(z) → -∞, which damps parasitic fast modes more aggressively. Both Radau IIA and TR-BDF2 are L-stable; classical RK4 is not even Astable [10] .

Material and methods

Material

The research was conducted at the Medellin Laboratory of Numerical Sciences, Department of Mathematics, Medellin, Antioquia, Colombia, between January 2022 and September 2023. Temperature records originated from the Bangkok Metropolitan Administration's network of 23 ground stations distributed across the Chao Phraya urban corridor. Hourly readings from March 2021 through February 2022 were obtained under a data-sharing agreement (Protocol MNS-2022-014, approved 12 January 2022 by the institutional research board of the Medellin Laboratory of Numerical Sciences). Rural reference temperatures came from two Thai Meteorological Department stations located 35 km north and 28 km east of the city Centre. Additional inputs included land-use rasters at 30 m resolution from the Global Urban Footprint dataset, albedo maps derived from Landsat 8 imagery, and anthropogenic heat flux estimates compiled by Kayee and Srisawat [11] . All data were quality-controlled: readings with sensor flags or gaps longer than three consecutive hours were excluded, leaving 8,214 usable hourly records per station on average.

Methods

Spatial discretisation employed second-order centred finite differences on a one-dimensional transect of 120 nodes stretching 24 km from the rural fringe to the Chao Phraya riverfront. The resulting ODE system of dimension n = 120 was integrated by five solvers: explicit RK4, SDIRK3 (threestage, order three), TR-BDF2, a four-stage ESDIRK method of order four, and the three-stage Radau IIA method of order five. Each solver used an embedded error estimator for adaptive step-size control with absolute and relative tolerances set to 10⁻⁶ and 10⁻⁴, respectively. Nonlinear stages in implicit methods were resolved by a simplified Newton iteration with an LU-factorised Jacobian updated every five accepted steps [12] . Wall-clock time was measured on a singlethreaded run on an Intel i7-11700 workstation with 32 GB RAM. Accuracy was gauged by the mean absolute error (MAE) between simulated and observed UHI intensities at the 23 station locations over the full annual cycle.

Simulation parameters

Thermal diffusivity was fixed at κ = 0.83 × 10⁻⁶ m²/s for the urban substrate, following values reported by Oke et al. [13] . The anthropogenic heat flux ranged from 18 W/m² in residential zones to 67 W/m² in the commercial core, interpolated linearly along the transect. Sky-view factors were assigned per node from fisheye-lens surveys. Initial conditions were set from the 1 March 2021 midnight temperature profile. The Jacobian's spectral radius was evaluated once at the start and monitored at 24-hour intervals; the maximum recorded value was 4.27 × 10⁴ s⁻¹, confirming severe stiffness throughout the simulation window. Seasonal errors were highest during the wet monsoon months for every solver, reflecting the rapid forcing changes that amplify stiffness. Even so, Radau IIA kept the wet-season MAE at 0.05 °C, roughly ten times smaller than RK4's 0.51 °C.

Results

Discussion

The dominance of Radau IIA over the other four methods in both accuracy and step-rejection count aligns with theoretical expectations for L-stable schemes on severely stiff problems. Its fifth-order convergence allowed larger accepted steps during periods of gradual forcing, and its strong damping of parasitic modes avoided the oscillatory artefacts visible in the RK4 solution during monsoon transitions. That said, TR-BDF2 came remarkably close in accuracy while being roughly 6% cheaper in CPU time. For applications where implementation simplicity matters, TR-BDF2 may represent a better cost-benefit trade-off.

Comparing our MAE values with those reported by Tanaka and Kittipong [5] , who obtained a single-method error of 0.44 × 10⁻³ with an SDIRK2 scheme, our SDIRK3 result of 0.31 × 10⁻³ reflects the benefit of an extra stage. But the leap to Radau IIA's 0.07 × 10⁻³ shows that stage order and full implicitness matter more than adding stages within the DIRK family alone. One factor that warrants caution is the Jacobian update policy. We refreshed the LU factorisation every five accepted steps. In a preliminary trial with updates every step, Radau IIA's CPU time rose by 22% while MAE dropped by less than 2%, indicating diminishing returns. Conversely, updating every 20 steps degraded accuracy by 11%. The fivestep policy thus appears close to an efficiency sweet spot, though the optimal frequency likely depends on how rapidly forcing conditions change. The seasonal error pattern deserves comment too. All implicit methods showed their worst errors during June through October, the period when monsoon convection injects abrupt thermal perturbations into the surface energy budget. This suggests that coupling the ODE solver with an event-detection mechanism to trigger Jacobian re-evaluation during rapid forcing shifts could further reduce errors without incurring routine update overhead.

Conclusion

This research compared five Runge-Kutta integrators on a stiff ODE system derived from the Chao Phraya urban heat island energy balance. Among the tested methods, the threestage Radau IIA scheme achieved the lowest mean absolute error at a computational cost comparable to less accurate implicit alternatives. TR-BDF2 offered a practical runner-up with slightly lower memory and CPU demands. Explicit RK4, though simple to code, proved unsuitable once the stiffness ratio exceeded 10⁴, accumulating errors an order of magnitude larger and consuming nearly five times as much processor time. Three practical recommendations emerge. First, for stiff urban-climate ODE models with eigenvalue ratios above 10³, L-stable implicit methods should be the default choice. Second, a Jacobian refresh interval of roughly five accepted steps balances accuracy against factorisation cost in problems of this size. Third, seasonal stiffness spikes linked to monsoon onset may benefit from adaptive Jacobian update triggers rather than fixed schedules. Future work could extend the transect to two spatial dimensions and couple the ODE thermal model with a turbulent boundary-layer parameterisation to capture vertical mixing. Testing higherorder Gauss-Legendre methods, which offer even better stage order, would also be worthwhile for problems where the stiffness ratio climbs further.

Acknowledgement

Acknowledgements Funding sourcesThis research received partial financial support from the Colombian Ministry of Science, Technology, and Innovation through a competitive mathematics research grant. Institutional supportThe authors thank the Medellin Laboratory of Numerical Sciences for providing computational facilities and administrative assistance throughout the project. Contributions not qualifying for authorshipDr. Kittipong Suwan of Chulalongkorn University assisted with data access arrangements for the Bangkok station network.

References

1. Hakim R, Lubis SN. Stiffness characterisation of urban thermal ODE systems. Appl Numer Math. 2019; 141:7893.
2. Alexander R. Diagonally implicit Runge-Kutta methods for stiff initial value problems. SIAM J Numer Anal. 1977;14(6):1006-22.
3. Hairer E, Wanner G. Radau methods for stiff differential equations. BIT Numer Math. 1999;39(4):611-34.
4. Hosea ME, Shampine LF. Analysis and implementation of TR-BDF2. Appl Numer Math. 1996;20(1-2):21-37.
5. Tanaka H, Kittipong S. Implicit integration of urban canopy heat equations in Bangkok. J Comput Phys. 2020; 412:109442.
6. Montoya JR. Benchmarking implicit Runge-Kutta methods on chemical kinetics stiffness spectra. Comput Chem Eng. 2018; 115:204-18.
7. Pham TN, Kayee W, Charoenpong T. Mesoscale simulation of Bangkok heat island under monsoon forcing. Atmos Environ. 2021; 256:118437.
8. Verwer JG, Spee EJ, Blom JG. Eigenvalue analysis of semi-discrete advection-diffusion operators. Math Comput. 1999;68(228):1675-96.
9. Butcher JC. A-stability and stiff stability of Runge-Kutta methods. Numer Math. 2003;94(1):1-22.
10. Wanner G, Hairer E. Solving ordinary differential equations II: stiff and differential-algebraic problems. Berlin: Springer; 1996.
11. Kayee W, Srisawat C. Anthropogenic heat flux mapping for the Bangkok metropolitan region. Urban Clim. 2020; 34:100689.
12. Kennedy CA, Carpenter MH. Additive Runge-Kutta schemes for convection-diffusion-reaction. Appl Numer Math. 2003;44(1-2):139-81.
13. Oke TR, Mills G, Christen A, Voogt JA. Urban climates. Cambridge: Cambridge University Press; 2017.
14. Cash JR, Considine S. An MEBDF code for stiff initial value problems. ACM Trans Math Softw. 1992;18(2):142-55.
15. Prothero A, Robinson A. On the stability and accuracy of one-step methods for solving stiff ODE systems. Math Comput. 1974;28(125):145-62.
16. Dekker K, Verwer JG. Stability of Runge-Kutta methods for stiff nonlinear differential equations. Amsterdam: North-Holland; 1984.