Skip to content

index.md

Latest commit

 

History

History
224 lines (187 loc) · 7.95 KB

index.md

File metadata and controls

224 lines (187 loc) · 7.95 KB

PositiveIntegrators.jl

The Julia library PositiveIntegrators.jl provides several time integration methods developed to preserve the positivity of numerical solutions.

Installation

PositiveIntegrators.jl is a registered Julia package. Thus, you can install it from the Julia REPL via

julia> using Pkg; Pkg.add("PositiveIntegrators")

If you want to update PositiveIntegrators.jl, you can use

julia> using Pkg; Pkg.update("PositiveIntegrators")

As usual, if you want to update PositiveIntegrators.jl and all other packages in your current project, you can execute

julia> using Pkg; Pkg.update()

Basic examples

Modified Patankar-Runge-Kutta schemes

Modified Patankar-Runge-Kutta (MPRK) schemes are unconditionally positive and conservative time integration schemes for the solution of positive and conservative ODE systems. The application of these methods is based on the representation of the ODE system as a so-called production-destruction system (PDS).

Production-destruction systems (PDS)

The application of MPRK schemes requires the ODE system to be represented as a production-destruction system (PDS). A PDS takes the general form

$$u_i'(t) = \sum_{j=1}^N \bigl(p_{ij}(t,\boldsymbol u) - d_{ij}(t,\boldsymbol u)\bigr),\quad i=1,\dots,N,$$

where \boldsymbol u=(u_1,\dots,u_n)^T is the vector of unknowns and both production terms p_{ij}(t,\boldsymbol u) and destruction terms d_{ij}(t,\boldsymbol u) must be nonnegative for all i,j=1,\dots,N. The meaning behind p_{ij} and d_{ij} is as follows:

  • p_{ij} with i\ne j represents the sum of all nonnegative terms which appear in equation i with a positive sign and in equation j with a negative sign.
  • d_{ij} with i\ne j represents the sum of all nonnegative terms which appear in equation i with a negative sign and in equation j with a positive sign.
  • p_{ii} represents the sum of all nonnegative terms which appear in equation i and don't have a negative counterpart in one of the other equations.
  • d_{ii} represents the sum of all negative terms which appear in equation i and don't have a positive counterpart in one of the other equations.

This naming convention leads to p_{ij} = d_{ji} for i≠ j and therefore a PDS is completely defined by the production matrix \mathbf{P}=(p_{ij})_{i,j=1,\dots,N} and the destruction vector \mathbf{d}=(d_{ii})_{i=1,\dots,N}.

As an example we consider the Lotka-Volterra model

$$\begin{aligned} u_1' &= 2u_1-u_1u_2,\\\ u_2' &= u_1u_2-u_2, \end{aligned}$$

which always has positive solutions if positive initial values are supplied. Assuming u_1,u_2>0, the above naming scheme results in

$$\begin{aligned} p_{11}(u_1,u_2) &= 2u_1,\\\ p_{21}(u_1,u_2) &= u_1u_2 = d_{12}(u_1,u_2) ,\\\ d_{22}(u_1,u_2) &= u_2, \end{aligned}$$

where all remaining production and destruction terms are zero. Consequently the production matrix \mathbf P and destruction vector \mathbf d are

$$\mathbf P(u_1,u_2) = \begin{pmatrix}2u_1 & 0\\ u_1u_2 & 0\end{pmatrix},\quad \mathbf d(u_1,u_2) = \begin{pmatrix}0\\ u_2\end{pmatrix}.$$ 
import Pkg; Pkg.add("OrdinaryDiffEq");  Pkg.add("Plots")

To solve this PDS together with initial values u_1(0)=u_2(0)=2 on the time domain (0,10), we first need to create a PDSProblem.

using PositiveIntegrators # load PDSProblem

function Pd(u, p, t)
  P = [2*u[1]  0; u[1]*u[2]  0] # Production matrix
  d = [0; u[2]] # Destruction vector
  return P, d
end

u0 = [2.0; 2.0] # initial values
tspan = (0.0, 10.0) # time span

# Create PDS
prob = PDSProblem(Pd, u0, tspan)
nothing #hide

Now that the problem has been created, we can solve it with any method of PositiveIntegrators.jl. In the following, we use the method MPRK22(1.0). In addition, we could also use any method provided by OrdinaryDiffEq.jl, but these might possibly generate negative approximations.

sol = solve(prob, MPRK22(1.0))
nothing # hide

Finally, we can use Plots.jl to visualize the solution.

using Plots

plot(sol)

Conservative production-destruction systems

A PDS with the additional property

$$p_{ii}(t,\boldsymbol y)=d_{ii}(t,\boldsymbol y)=0$$

for i=1,\dots,N is called conservative. In this case we have p_{ij}=d_{ji} for all i,j=1,\dots,N, which leads to

$$\frac{d}{dt}\sum_{i=1}^N y_i=\sum_{i=1}^N y_i' = \sum_{\mathclap{i,j=1}}^N \bigl(p_{ij}(t,\boldsymbol y) - d_{ij}(t,\boldsymbol y)\bigr)= \sum_{\mathclap{i,j=1}}^N \bigl(p_{ij}(t,\boldsymbol y) - p_{ji}(t,\boldsymbol y)\bigr) = 0.$$

This shows that the sum of the state variables of a conservative PDS remains constant over time, i.e.

$$\sum_{i=1}^N y_i(t) = \sum_{i=1}^N y_i(0)$$

for all times t>0.

One specific example of a conservative PDS is the SIR model

$$S' = -\frac{β S I}{N},\quad I'= \frac{β S I}{N} - γ I,\quad R'=γ I,$$

with N=S+I+R and \beta,\gamma>0. Assuming S,I,R>0 the production and destruction terms are given by

$$p_{21}(S,I,R) = d_{12}(S,I,R) = \frac{β S I}{N},\quad p_{32}(S,I,R) = d_{23}(S,I,R) = γ I,$$

where the remaining production and destruction terms are zero. The corresponding production matrix \mathbf P is

$$\mathbf P(S,I,R) = \begin{pmatrix}0 & 0 & 0\\ \frac{β S I}{N} & 0 & 0\\ 0 & γ I & 0\end{pmatrix}.$$

The following example shows how to implement the above SIR model with \beta=0.4, \gamma=0.04, initial conditions S(0)=997, I(0)=3, R(0)=0 and time domain (0, 100) using ConservativePDSProblem from PositiveIntegrators.jl.

import Pkg; Pkg.add("OrdinaryDiffEq");

using PositiveIntegrators

# Out-of-place implementation of the P matrix for the SIR model
function P(u, p, t)
  S, I, R = u

  β = 0.4
  γ = 0.04
  N = 1000.0

  P = zeros(3,3)
  P[2,1] = β*S*I/N
  P[3,2] = γ*I
  return P
end

u0 = [997.0; 3.0; 0.0]; # initial values
tspan = (0.0, 100.0); # time span

# Create SIR problem
prob = ConservativePDSProblem(P, u0, tspan)
nothing # hide

Since the SIR model is not only conservative but also positive, we can use any scheme from PositiveIntegrators.jl to solve it. Here we use MPRK22(1.0). Please note that any method from OrdinaryDiffEq.jl can be used as well, but might possibly generate negative approximations.

sol = solve(prob, MPRK22(1.0))
nothing # hide

Finally, we can use Plots.jl to visualize the solution.

using Plots

plot(sol, label = ["S" "I" "R"], legend=:right)
plot!(sol.t, sum.(sol.u), label = "S+I+R") # Plot S+I+R over time.

We see that there is always a nonnegative number of people in each compartment, while the population S+I+R remains constant over time.

Referencing

If you use PositiveIntegrators.jl for your research, please cite it using the bibtex entry

@misc{PositiveIntegrators.jl,
  title={{PositiveIntegrators.jl}: {A} {J}ulia library of positivity-preserving
         time integration methods},
  author={Kopecz, Stefan and Ranocha, Hendrik and contributors},
  year={2023},
  doi={10.5281/zenodo.10868393},
  url={https://github.com/SKopecz/PositiveIntegrators.jl}
}

License and contributing

This project is licensed under the MIT license (see License). Since it is an open-source project, we are very happy to accept contributions from the community. Please refer to the section Contributing for more details.