GitHub

Tutorial: Solution of Robertson problem

In this tutorial we show that MPRK schemes can be used to integrate stiff problems. We also show how callbacks can be used to change the time step in non-adaptive schemes.

Definition of the production-destruction system

The well known Robertson problem is given by

\[\begin{aligned} u_1' &= -0.04u_1+10^4 u_2u_3, & u_1(0)&=1,\\ u_2' &= 0.04u_1-10^4 u_2u_3-3⋅10^7 u_2^2, & u_2(0)&=0 \\ u_3' &= 3⋅10^7 u_2^2, & u_3(0)&=0. \end{aligned}\]

The time domain of interest is $t\in[0,10^{11}]$, because of which some kind of adaptive time stepping is required.

The model can be represented as a conservative PDS with production terms

\[\begin{aligned} p_{12}(t,\mathbf{u}) &= 10^4u_2u_3,& p_{21}(t,\mathbf{u}) &= 0.04u_1, & p_{32}(t,\mathbf{u}) &= 3⋅10^7u_2^2, \end{aligned}\]

whereby production terms not listed have the value zero. Since the PDS is conservative, we have $d_{ij}=p_{ji}$ and the system is fully determined by the production matrix $\mathbf P=(p_{ij})$.

Solution of the production-destruction system

Now we are ready to define a ConservativePDSProblem and to solve this problem with any method of PositiveIntegrators.jl or OrdinaryDiffEq.jl which is suited for stiff problems.

Since this PDS consists of only three differential equations we provide an out-of-place implementation for the production matrix. Furthermore, we use static arrays from StaticArrays.jl for additional efficiency. See also the tutorials on the solution of an NPZD model or an stratospheric reaction problem.

using PositiveIntegrators, StaticArraysfunction prod(u, p, t)    @SMatrix [0.0 1e4*u[2]*u[3] 0.0;              4e-2*u[1] 0.0 0.0;              0.0 3e7*u[2]^2 0.0]endu0 = @SVector [1.0, 0.0, 0.0]tspan = (0.0, 1.0e11)prob = ConservativePDSProblem(prod, u0, tspan)sol = solve(prob, MPRK43I(1.0, 0.5))
using Plotsplot(sol, tspan = (1e-6, 1e11),  xaxis = :log,     idxs = [(0, 1), ((x, y) -> (x, 1e4 .* y), 0, 2), (0, 3)],     label = ["u₁" "10⁴u₂" "u₃"])
Example block output

PositiveIntegrators.jl provides the function isnonnegative (and also isnegative) to check if the solution is actually nonnegative, as expected from an MPRK scheme.

isnonnegative(sol)
true

Using callbacks to solve the Robertson problem with non-adaptive schemes

The SSPMPRK43() scheme is only available with fixed time stepping. With a scheme like this, it would take a huge amount of time to solve the Robertson problem, since the time step must be chosen very small to accurately solve the problem in its initial phase. However, the use of a callback allows us to modify the time step size after each step, which makes a solution with a fixed step method possible.

In the following example the callback increases the time step size by a factor of 1.5 after each time step.

using DiffEqCallbacksusing DiffEqBasestepsize_callback = DiscreteCallback(    Returns(true), # adapt the step size after every time step    integrator -> set_proposed_dt!(integrator, 1.5 * get_proposed_dt(integrator));    save_positions = (false, false),    initialize = (c, u, t, integrator) -> set_proposed_dt!(integrator, 1.0e-5))sol_cb = solve(prob, SSPMPRK43(); dt = Inf, callback =  stepsize_callback);
plot(sol_cb, tspan = (1e-6, 1e11),  xaxis = :log,     idxs = [(0, 1), ((x, y) -> (x, 1e4 .* y), 0, 2), (0, 3)],     label = ["u₁" "10⁴u₂" "u₃"])
Example block output

This solution is also nonnegative.

isnonnegative(sol_cb)
true

Package versions

These results were obtained using the following versions.

using InteractiveUtilsversioninfo()println()using PkgPkg.status(["PositiveIntegrators", "StaticArrays", "LinearSolve", "DiffEqCallbacks", "DiffEqBase"],           mode=PKGMODE_MANIFEST)
Julia Version 1.13.0
Commit d1c37793dd2 (2026-09-09 19:00 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 9V74 80-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)
Environment:
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager

Status `~/work/PositiveIntegrators.jl/PositiveIntegrators.jl/docs/Manifest.toml`
  [14f7f29c] AMD v0.5.4
  [4fba245c] ArrayInterface v7.30.2
  [70df07ce] BracketingNonlinearSolve v1.12.7
  [2569d6c7] ConcreteStructs v0.2.8
  [864edb3b] DataStructures v0.19.6
  [2b5f629d] DiffEqBase v7.21.2
  [459566f4] DiffEqCallbacks v4.19.4
  [a0c0ee7d] DifferentiationInterface v0.7.21
  [ffbed154] DocStringExtensions v0.9.5
  [4e289a0a] EnumX v1.0.7
  [7034ab61] FastBroadcast v1.4.0
  [9aa1b823] FastClosures v0.3.2
  [a4df4552] FastPower v1.5.0
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v1.13.0
  [46192b85] GPUArraysCore v0.2.0
  [ba0b0d4f] Krylov v0.10.10
  [2faa5264] LHLFactorization v2.2.2
  [7ed4a6bd] LinearSolve v5.18.0
  [46d2c3a1] MuladdMacro v0.2.7
  [bbf590c4] OrdinaryDiffEqCore v4.18.0
  [d1b20bf0] PositiveIntegrators v0.2.21 `~/work/PositiveIntegrators.jl/PositiveIntegrators.jl`
  [d236fae5] PreallocationTools v1.7.1
  [aea7be01] PrecompileTools v1.3.4
  [21216c6a] Preferences v1.6.0
  [0c0d3e7f] PureKLU v1.6.0
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v4.5.1
  [189a3867] Reexport v1.2.2
  [9fe22ead] RespecializeParams v1.3.0
  [0bca4576] SciMLBase v3.55.0
  [a6db7da4] SciMLLogging v2.1.0
  [c0aeaf25] SciMLOperators v1.30.1
  [53ae85a6] SciMLStructures v1.10.5
  [efcf1570] Setfield v1.1.2
  [a57abbd0] SparseColumnPivotedQR v2.1.8
  [90137ffa] StaticArrays v1.9.22
  [1e83bf80] StaticArraysCore v1.4.4
  [10745b16] Statistics v1.11.5
  [2efcf032] SymbolicIndexingInterface v0.3.55
  [781d530d] TruncatedStacktraces v1.4.0
  [856f044c] MKL_jll v2025.2.0+0
  [b77e0a4c] InteractiveUtils v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.13.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [de0858da] Printf v1.11.0
  [9a3f8284] Random v1.11.0
  [2f01184e] SparseArrays v1.13.0
  [4536629a] OpenBLAS_jll v0.3.30+0