Tutorial: Solving a scalar production-destruction equation with MPRK schemes

Originally, modified Patankar-Runge-Kutta (MPRK) schemes were designed to solve positive and conservative systems of ordinary differential equations. The conservation property requires that the system consists of at least two scalar differential equations. Nevertheless, we can also apply the idea of the Patankar trick to a scalar production-destruction system (PDS)

\[u'(t)=p(u(t))-d(u(t)),\quad u(0)=u_0>0\]

with nonnegative functions $p$ and $d$. Since conservation is not an issue here, we can apply the Patankar trick to the destruction term $d$ to ensure positivity and leave the production term $p$ unweighted. A first-order scheme of this type, based on the forward Euler method, reads

\[u^{n+1}= u^n + Δ t p(u^n) - Δ t d(u^n)\frac{u^{n+1}}{u^n}\]

and this idea can easily be generalized to higher-order explicit Runge-Kutta schemes.

By closer inspection we realize that this is exactly the approach the MPRK schemes of PositiveIntegrators.jl use to solve non-conservative PDS for which the production matrix is diagonal. Hence, we can use the existing schemes to solve a scalar PDS by regarding the production term as a $1×1$-matrix and the destruction term as a $1$-vector.

Example 1

We want to solve

\[u' = u^2 - u,\quad u(0) = 0.95\]

for $0≤ t≤ 10$. Here, we can choose $p(u)=u^2$ as production term and $d(u)=u$ as destruction term. The exact solution of this problem is

\[u(t) = \frac{19}{19+e^{t}}.\]

Next, we show how to solve this scalar PDS in the way discussed above. Please note that we must use PDSProblem to create the problem. Furthermore, we use static matrices and vectors from StaticArrays.jl instead of standard arrays for efficiency.

using PositiveIntegrators, StaticArrays, Plots

u0 = @SVector [0.95] # 1-vector
tspan = (0.0, 10.0)

# Attention: Input u is a 1-vector
prod(u, p, t) = @SMatrix [u[1]^2] # create static 1x1-matrix
dest(u, p, t) = @SVector [u[1]] # create static 1-vector
prob = PDSProblem(prod, dest, u0, tspan)

sol = solve(prob, MPRK22(1.0))

# plot
tt = 0:0.1:10
f(t) = 19.0 / (19.0 + exp(t)) # exact solution
plot(tt, f.(tt), label="exact")
plot!(sol, label="u")

Example 2

Next, we want to compute positive solutions of a more challenging scalar PDS. In Example 1, we could have also used standard schemes from OrdinaryDiffEq.jl and use the solver option isoutofdomain to ensure positivity. But this is not always the case as the following example will show.

We want to compute the nonnegative solution of

\[u'(t) = -\sqrt{\lvert u(t)\rvert },\quad u(0)=1.\]

for $t≥ 0$. Please note that this initial value problem has infinitely many solutions

\[u(t) = \begin{cases} \frac{1}{4}(t-2)^2, & 0≤ t< 2,\\ 0, & 2≤ t < t^*,\\ -\frac{1}{4}(t-2)^2, & t^*≤ t, \end{cases}\]

where $t^*≥ 2$ is arbitrary. But among these, the only nonnegative solution is

\[u(t) = \begin{cases} \frac{1}{4}(t-2)^2, & 0≤ t< 2,\\ 0, & 2≤ t. \end{cases}\]

This is the solution we want to compute.

First, we try this using a standard solver from OrdinaryDiffEq.jl. We try to enforce positivity with the solver option isoutofdomain by specifying that negative solution components are not acceptable.

using OrdinaryDiffEqRosenbrock

tspan = (0.0, 3.0)
u0 = 1.0

f(u, p, t) = -sqrt(abs(u))
prob = ODEProblem(f, u0, tspan)

sol = solve(prob, Rosenbrock23(); isoutofdomain = (u, p, t) -> any(<(0), u))

retcode: Unstable
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 83-element Vector{Float64}:
 0.0
 0.003163858403911277
 0.034802442443024044
 0.09367771717181361
 0.17538247876115193
 0.289520297678292
 0.4225451283138154
 0.5773542035170329
 0.7383044438936105
 0.904052247269852
 ⋮
 1.9960526866394865
 1.9960526866395492
 1.9960526866395742
 1.9960526866395842
 1.9960526866395882
 1.9960526866395962
 1.9960526866396122
 1.9960526866396135
 1.996052686639614
u: 83-element Vector{Float64}:
 1.0
 0.9968386434463932
 0.965499712213253
 0.908511318269506
 0.8322910840841042
 0.7313873654216926
 0.621992024565334
 0.5057980826645783
 0.3976922095042322
 0.2998942510751979
 ⋮
 3.8149886585940746e-27
 8.646781239421471e-28
 2.763213158726151e-28
 1.3389541665994338e-28
 9.146973148338877e-29
 2.999872984387982e-29
 2.1464764270874556e-30
 6.5547154986337885e-31
 3.0324240953021175e-31

We see that isoutofdomain cannot be used to ensure nonnegative solutions in this case, as the computation stops at about $t≈ 2$ before the desired final time is reached. For at least first- and second-order explicit Runge-Kutta schemes, this can also be shown analytically. A brief computation reveals that to ensure nonnegative solutions, the time step size must tend to zero if the numerical solution tends to zero.

Next, we want to use an MPRK scheme. We can choose $p(u)=0$ as the production term and $d(u)=\sqrt{\lvert u\rvert }$ as the destruction term. Furthermore, we create the PDSProblem in the same way as in Example 1.

using PositiveIntegrators, StaticArrays, Plots

tspan = (0.0, 3.0)
u0 = @SVector [1.0]

prod(u, p, t) = @SMatrix zeros(1,1)
dest(u, p, t) = @SVector [sqrt(abs(first(u)))]
prob = PDSProblem(prod, dest, u0, tspan)

sol = solve(prob, MPRK22(1.0))

# plot
tt = 0:0.03:3
f(t) = 0.25 * (t - 2)^2 * (t <= 2) # exact solution
plot(tt, f.(tt), label="exact")
plot!(sol, label="u")

We can see that the MPRK scheme used is well suited to solve the problem.

Package versions

These results were obtained using the following versions.

using InteractiveUtils
versioninfo()
println()

using Pkg
Pkg.status(["PositiveIntegrators", "StaticArrays", "OrdinaryDiffEqRosenbrock"],
           mode = PKGMODE_MANIFEST)

Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 7763 64-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager

Status `~/work/PositiveIntegrators.jl/PositiveIntegrators.jl/docs/Manifest.toml`
  [43230ef6] OrdinaryDiffEqRosenbrock v1.10.0
  [d1b20bf0] PositiveIntegrators v0.2.12 `~/work/PositiveIntegrators.jl/PositiveIntegrators.jl`
  [90137ffa] StaticArrays v1.9.13