Tutorial: Positive-projection method
This tutorial is about solving an ODE using the projection method introduced by Adrian Sandu in Positive Numerical Integration Methods for Chemical Kinetic Systems. It guarantees positivity by solving an optimization problem while preserving all linear invariants.
The Sandu projection is a post-processing technique that can be used in combination with any ODE solver. If the ODE solver computes a negative approximation at any time step, the projection method calculates a positive approximation, also taking into account the linear invariants.
Solution of the ODE system
As an example we want to the solve the NPZD problem prob_pds_npzd, which is an ODE system in which negative approximations quickly lead to unacceptable solutions. First, we solve the problem without Sandu projection and select ROS2 form OrdinaryDiffEq.jl as ODE solver.
using PositiveIntegratorsusing OrdinaryDiffEqRosenbrockusing Plotsprob = prob_pds_npzdref_sol = solve(prob, ROS2(); abstol = 1e-8, reltol = 1e-6); # reference solution for plottingsol = solve(prob, ROS2(); abstol = 5e-2, reltol = 1e-1)plot(ref_sol, linestyle = :dash, label = "", color = palette(:default)[1:4]', plotdensity = 1000)plot!(sol, ylims = (-2.5, 12.5), denseplot = false, markers = :circle, linewidth = 2, color = palette(:default)[1:4]', label = ["N" "P" "Z" "D"], legend = :right)
The plot shows the numerical solution obtained with ROS2 compared to a reference solution (dashed lines). We see that the ROS2 method produces negative approximations, which can occur because Rosenbrock methods are not positivity-preserving. For the NPZD problem, however, this is fatal and leads to a completely unacceptable numerical solution. It is therefore particularly important to use techniques that guarantee positivity of the numerical approximations for this problem. We achieve this below with the SanduProjection.
To apply the SanduProjection we need to choose an optimization solver which is supported by JuMP.jl and can handle quadratic optimization problems (QP). In this tutorial we select Clarabel.jl as optimization solver.
In addition, we need to specify the linear invariants of the problem. The only linear invariant of the NPZD problem is $N(t)+P(t)+Z(t)+D(t)=N(0)+P(0)+Z(0)+D(0)=15$ for all times $t≥0$. This can be written in the form
\[\mathbf{A}^T \begin{pmatrix} N(t)\\ P(t)\\ Z(t)\\ D(t) \end{pmatrix} = \mathbf{b}\]
with $\mathbf{A}^T = [1.0,\ 1.0,\ 1.0,\ 1.0]$ and $\mathbf{b} = [15]$.
The projection method SanduProjection is implemented as a callback and hence, must be passed as an argument to the keyword callback. In addition, we must also use save_everystep = false.
using JuMP, ClarabelAT = [1.0 1.0 1.0 1.0]b = [15.0]proj = SanduProjection(Model(Clarabel.Optimizer), AT, b)sol_proj = solve(prob, ROS2(); abstol = 5e-2, reltol = 1e-1, save_everystep = false, callback = proj);plot(ref_sol, linestyle = :dash, label = "", color = palette(:default)[1:4]', plotdensity = 1000)plot!(sol_proj, ylims = (-2.5, 12.5), denseplot = false, markers = :circle, linewidth = 2, color = palette(:default)[1:4]', label = ["N" "P" "Z" "D"], legend = :right)
As intended, negative approximations no longer occur and we obtain an acceptable approximation.
The SanduProjection is implemented as a DiscreteCallback and we can display the number of projection steps in the following way.
@show get_numsteps_SanduProjection(proj)1We can see that in this example, a single projection step was already sufficient.
Package versions
These results were obtained using the following versions.
using InteractiveUtilsversioninfo()println()using PkgPkg.status(["PositiveIntegrators", "JuMP", "Clarabel", "OrdinaryDiffEqRosenbrock", "Plots"], 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`
[47edcb42] ADTypes v1.24.0
[14f7f29c] AMD v0.5.4
[4fba245c] ArrayInterface v7.30.2
[61c947e1] Clarabel v0.11.1
[d38c429a] Contour v0.6.3
[864edb3b] DataStructures v0.19.6
[2b5f629d] DiffEqBase v7.21.2
[a0c0ee7d] DifferentiationInterface v0.7.21
[c87230d0] FFMPEG v0.4.5
[7034ab61] FastBroadcast v1.4.0
[6a86dc24] FiniteDiff v2.33.0
⌅ [53c48c17] FixedPointNumbers v0.8.6
[f6369f11] ForwardDiff v1.4.6
[28b8d3ca] GR v0.73.27
⌅ [14197337] GenericLinearAlgebra v0.3.19
[1019f520] JLFzf v0.1.11
[682c06a0] JSON v1.9.0
[4076af6c] JuMP v1.31.2
[b964fa9f] LaTeXStrings v1.4.1
[23fbe1c1] Latexify v0.16.12
[7ed4a6bd] LinearSolve v5.18.0
[1914dd2f] MacroTools v0.5.16
[b8f27783] MathOptInterface v1.54.0
[442fdcdd] Measures v0.3.3
[46d2c3a1] MuladdMacro v0.2.7
[d8a4904e] MutableArithmetics v1.8.1
[77ba4419] NaNMath v1.1.4
[bac558e1] OrderedCollections v2.0.1
[bbf590c4] OrdinaryDiffEqCore v4.18.0
[4302a76b] OrdinaryDiffEqDifferentiation v3.12.3
[43230ef6] OrdinaryDiffEqRosenbrock v2.7.4
[b4bd8bb3] OrdinaryDiffEqRosenbrockTableaus v2.4.2
[ccf2f8ad] PlotThemes v3.3.0
[995b91a9] PlotUtils v1.5.0
[91a5bcdd] Plots v1.41.7
[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
[bfc457fd] QDLDL v0.4.1
[3cdcf5f2] RecipesBase v1.3.4
[01d81517] RecipesPipeline v0.6.12
[731186ca] RecursiveArrayTools v4.5.1
[189a3867] Reexport v1.2.2
[05181044] RelocatableFolders v1.0.1
[ae029012] Requires v1.3.1
[0bca4576] SciMLBase v3.55.0
[6c6a2e73] Scratch v1.3.0
[992d4aef] Showoff v1.1.1
[66db9d55] SnoopPrecompile v1.0.3
[90137ffa] StaticArrays v1.9.22
[10745b16] Statistics v1.11.5
[2913bbd2] StatsBase v0.34.13
[2efcf032] SymbolicIndexingInterface v0.3.55
⌅ [a759f4b9] TimerOutputs v0.5.29
[1cfade01] UnicodeFun v0.4.1
[41fe7b60] Unzip v0.2.0
[2a0f44e3] Base64 v1.11.0
[ade2ca70] Dates v1.11.0
[f43a241f] Downloads v1.7.0
[37e2e46d] LinearAlgebra v1.13.0
[44cfe95a] Pkg v1.13.0
[de0858da] Printf v1.11.0
[3fa0cd96] REPL v1.11.0
[9a3f8284] Random v1.11.0
[2f01184e] SparseArrays v1.13.0
[4607b0f0] SuiteSparse
[fa267f1f] TOML v1.0.3
[cf7118a7] UUIDs v1.11.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`