PDSProblemLibrary.jl

# linear model problem

P_linmod(u, p, t) = @SMatrix [0.0 u[2]; 5.0*u[1] 0.0]
f_linmod(u, p, t) = @SVector [u[2] - 5.0 * u[1]; 5.0 * u[1] - u[2]]
function f_linmod_analytic(u0, p, t)
    u₁⁰, u₂⁰ = u0
    a = 5.0
    b = 1.0
    c = a + b
    return ((u₁⁰ + u₂⁰) * [b; a] + exp(-c * t) * (a * u₁⁰ - b * u₂⁰) * [1; -1]) / c
end
u0_linmod = @SVector [0.9, 0.1]
"""
    prob_pds_linmod

Positive and conservative autonomous linear PDS
```math
\\begin{aligned}
u_1' &= u_2 - 5u_1,\\\\
u_2' &= 5u_1 - u_2,
\\end{aligned}
```
with initial value ``\\mathbf{u}_0 = (0.9, 0.1)^T`` and time domain ``(0.0, 2.0)``.
There is one independent linear invariant, e.g. ``u_1+u_2 = 1``.

## References

- Hans Burchard, Eric Deleersnijder, and Andreas Meister.
  "A high-order conservative Patankar-type discretisation for stiff systems of
  production-destruction equations."
  Applied Numerical Mathematics 47.1 (2003): 1-30.
  [DOI: 10.1016/S0168-9274(03)00101-6](https://doi.org/10.1016/S0168-9274(03)00101-6)
"""
prob_pds_linmod = ConservativePDSProblem(P_linmod, u0_linmod, (0.0, 2.0),
                                         analytic = f_linmod_analytic, std_rhs = f_linmod,
                                         linear_invariants = @SMatrix[1.0 1.0])

function P_linmod!(P, u, p, t)
    P .= P_linmod(u, p, t)
    return nothing
end
function f_linmod!(du, u, p, t)
    du .= f_linmod(u, p, t)
    return nothing
end

"""
    prob_pds_linmod_inplace

Same as [`prob_pds_linmod`](@ref) but with in-place computation.
"""
prob_pds_linmod_inplace = ConservativePDSProblem(P_linmod!, Array(u0_linmod),
                                                 (0.0, 2.0),
                                                 analytic = f_linmod_analytic,
                                                 std_rhs = f_linmod!,
                                                 linear_invariants = @SMatrix[1.0 1.0])

# nonlinear model problem
function P_nonlinmod(u, p, t)
    @SMatrix [0.0 0.0 0.0; u[2] * u[1]/(u[1] + 1.0) 0.0 0.0; 0.0 0.3*u[2] 0.0]
end
function f_nonlinmod(u, p, t)
    @SVector [-u[2] * u[1] / (u[1] + 1);
              u[2] * u[1] / (u[1] + 1) - 0.3 * u[2];
              0.3 * u[2]]
end
u0_nonlinmod = @SVector [9.98; 0.01; 0.01]
"""
    prob_pds_nonlinmod

Positive and conservative autonomous nonlinear PDS
```math
\\begin{aligned}
u_1' &= -\\frac{u_1u_2}{u_1 + 1.0},\\\\
u_2' &= \\frac{u_1u_2}{u_1 + 1.0} - 0.3u_2,\\\\
u_3' &= 0.3 u_2,
\\end{aligned}
```
with initial value ``\\mathbf{u}_0 = (9.98, 0.01, 0.01)^T`` and time domain ``(0.0, 30.0)``.
There is one independent linear invariant, e.g. ``u_1+u_2+u_3 = 10.0``.

## References

- Hans Burchard, Eric Deleersnijder, and Andreas Meister.
  "A high-order conservative Patankar-type discretisation for stiff systems of
  production-destruction equations."
  Applied Numerical Mathematics 47.1 (2003): 1-30.
  [DOI: 10.1016/S0168-9274(03)00101-6](https://doi.org/10.1016/S0168-9274(03)00101-6)
"""
prob_pds_nonlinmod = ConservativePDSProblem(P_nonlinmod, u0_nonlinmod, (0.0, 30.0),
                                            std_rhs = f_nonlinmod,
                                            linear_invariants = @SMatrix[1.0 1.0 1.0])

# robertson problem
function P_robertson(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]
end
function f_robertson(u, p, t)
    return @SVector [1e4 * u[2] * u[3] - 4e-2 * u[1];
                     4e-2 * u[1] - 1e4 * u[2] * u[3] - 3e7 * u[2]^2;
                     3e7 * u[2]^2]
end;
u0_robertson = @SVector [1.0, 0.0, 0.0]
"""
    prob_pds_robertson

Positive and conservative autonomous nonlinear PDS
```math
\\begin{aligned}
u_1' &= -0.04u_1+10^4 u_2u_3,\\\\
u_2' &=  0.04u_1-10^4 u_2u_3-3⋅10^7 u_2^2,\\\\
u_3' &= 3⋅10^7 u_2^2,
\\end{aligned}
```
with initial value ``\\mathbf{u}_0 = (1.0, 0.0, 0.0)^T`` and time domain ``(0.0, 10^{11})``.
There is one independent linear invariant, e.g. ``u_1+u_2+u_3 = 1.0``.

## References

- Ernst Hairer, Gerd Wanner.
  "Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems."
  2nd Edition, Springer (2002): Section IV.1.
"""
prob_pds_robertson = ConservativePDSProblem(P_robertson, u0_robertson, (0.0, 1.0e11),
                                            std_rhs = f_robertson,
                                            linear_invariants = @SMatrix[1.0 1.0 1.0])

# brusselator problem
function P_brusselator(u, p, t)
    u2u5 = u[2] * u[5]

    @SMatrix [0.0 0.0 0.0 0.0 0.0 0.0;
              0.0 0.0 0.0 0.0 0.0 0.0;
              0.0 u2u5 0.0 0.0 0.0 0.0;
              0.0 0.0 0.0 0.0 u[5] 0.0;
              u[1] 0.0 0.0 0.0 0.0 u[5]^2*u[6];
              0.0 0.0 0.0 0.0 u2u5 0.0]
end;
function f_brusselator(u, p, t)
    @SVector [-u[1];
              -u[2] * u[5];
              u[2] * u[5];
              u[5];
              u[1] - u[2] * u[5] + u[5]^2 * u[6] - u[5];
              u[2] * u[5] - u[5]^2 * u[6]]
end;
u0_brusselator = @SVector [10.0, 10.0, 0.0, 0.0, 0.1, 0.1]
"""
    prob_pds_brusselator

Positive and conservative autonomous nonlinear PDS
```math
\\begin{aligned}
u_1' &= -u_1,\\\\
u_2' &= -u_2u_5,\\\\
u_3' &= u_2u_5,\\\\
u_4' &= u_5,\\\\
u_5' &= u_1 - u_2u_5 + u_5^2u_6 - u_5,\\\\
u_6' &= u_2u_5 - u_5^2u_6,
\\end{aligned}
```
with initial value ``\\mathbf{u}_0 = (10.0, 10.0, 0.0, 0.0, 0.1, 0.1)^T`` and time domain ``(0.0, 20.0)``.
There are two independent linear invariants, e.g. ``u_1+u_4+u_5+u_6 = 10.2`` and ``u_2+u_3 = 10.0``.

## References

- Luca Bonaventura,  and Alessandro Della Rocca.
  "Unconditionally Strong Stability Preserving Extensions of the TR-BDF2 Method."
  Journal of Scientific Computing 70 (2017): 859 - 895.
  [DOI: 10.1007/s10915-016-0267-9](https://doi.org/10.1007/s10915-016-0267-9)
"""
prob_pds_brusselator = ConservativePDSProblem(P_brusselator, u0_brusselator, (0.0, 10.0),
                                              std_rhs = f_brusselator,
                                              linear_invariants = @SMatrix[1.0 0.0 0.0 1.0 1.0 1.0;
                                                                           0.0 1.0 1.0 0.0 0.0 0.0])

# SIR problem
P_sir(u, p, t) = @SMatrix [0.0 0.0 0.0; 2*u[1]*u[2] 0.0 0.0; 0.0 u[2] 0.0]
f_sir(u, p, t) = @SVector [-2 * u[1] * u[2];
                           2 * u[1] * u[2] - u[2];
                           u[2]];
u0_sir = @SVector [0.99, 0.005, 0.005]
"""
    prob_pds_sir

Positive and conservative autonomous nonlinear PDS
```math
\\begin{aligned}
u_1' &= -2u_1u_2,\\\\
u_2' &= 2u_1u_2 - u_2,\\\\
u_3' &= u_2,
\\end{aligned}
```
with initial value ``\\mathbf{u}_0 = (0.99, 0.005, 0.005)^T`` and time domain ``(0.0, 20.0)``.
There is one independent linear invariant, e.g. ``u_1+u_2+u_3 = 1.0``.

## References

- Ronald E. Mickens, and Talitha M. Washington.
  "NSFD discretizations of interacting population models satisfying conservation laws."
  Computers and Mathematics with Applications 66 (2013): 2307-2316.
  [DOI: 10.1016/j.camwa.2013.06.011](https://doi.org/10.1016/j.camwa.2013.06.011)
"""
prob_pds_sir = ConservativePDSProblem(P_sir, u0_sir, (0.0, 20.0), std_rhs = f_sir,
                                      linear_invariants = @SMatrix[1.0 1.0 1.0])

# bertolazzi problem
function P_bertolazzi(u, p, t)
    f1 = 5 * u[2] * u[3] / (1e-2 + (u[2] * u[3])^2) +
         u[2] * u[3] / (1e-16 + u[2] * u[3] * (1e-8 + u[2] * u[3]))
    f2 = 10 * u[1] * u[3]^2
    f3 = 0.1 * (u[3] - u[2] - 2.5)^2 * u[1] * u[2]

    return @SMatrix [0.0 f1 f1; f2 0.0 f2; f3 f3 0.0]
end
function f_bertolazzi(u, p, t)
    f1 = 5 * u[2] * u[3] / (1e-2 + (u[2] * u[3])^2) +
         u[2] * u[3] / (1e-16 + u[2] * u[3] * (1e-8 + u[2] * u[3]))
    f2 = 10 * u[1] * u[3]^2
    f3 = 0.1 * (u[3] - u[2] - 2.5)^2 * u[1] * u[2]

    return @SVector [2 * f1 - f2 - f3; -f1 + 2 * f2 - f3; -f1 - f2 + 2 * f3]
end
u0_bertolazzi = @SVector [0.0, 1.0, 2.0]
"""
    prob_pds_bertolazzi

Positive and conservative autonomous nonlinear PDS
```math
\\begin{aligned}
\\mathbf{u}'=\\begin{pmatrix}2 &-1 &-1\\\\-1 &2 &-1\\\\-1& -1& 2\\end{pmatrix}\\begin{pmatrix}5u_2u_3/(10^{-2} + (u_2u_3)^2) + u_2u_3/(10^{-16} + u_2u_3(10^{-8} + u_2u_3))\\\\
10u_1u_3^2\\\\
0.1(u_3 - u_2 - 2.5)^2u_1u_2\\end{pmatrix}
\\end{aligned}
```
with initial value ``\\mathbf{u}_0 = (0.0, 1.0, 2.0)^T`` and time domain ``(0.0, 1.0)``.
There is one independent linear invariant, e.g. ``u_1+u_2+u_3 = 3.0``.

## References

- Enrico Bertolazzi.
  "Positive and conservative schemes for mass action kinetics."
  Computers and Mathematics with Applications 32 (1996): 29-43.
  [DOI: 10.1016/0898-1221(96)00142-3](https://doi.org/10.1016/0898-1221(96)00142-3)
"""
prob_pds_bertolazzi = ConservativePDSProblem(P_bertolazzi, u0_bertolazzi, (0.0, 1.0),
                                             std_rhs = f_bertolazzi,
                                             linear_invariants = @SMatrix[1.0 1.0 1.0])

# npzd problem
function P_npzd(u, p, t)
    dnp = u[1] / (0.01 + u[1]) * u[2]
    dpz = 0.5 * (1.0 - exp(-1.21 * u[2]^2)) * u[3]
    dpn = 0.01 * u[2]
    dzn = 0.01 * u[3]
    ddn = 0.003 * u[4]
    dpd = 0.05 * u[2]
    dzd = 0.02 * u[3]

    return @SMatrix [0.0 dpn dzn ddn; dnp 0.0 0.0 0.0; 0.0 dpz 0.0 0.0; 0.0 dpd dzd 0.0]
end
function f_npzd(u, p, t)
    dnp = u[1] / (0.01 + u[1]) * u[2]
    dpz = 0.5 * (1.0 - exp(-1.21 * u[2]^2)) * u[3]
    dpn = 0.01 * u[2]
    dzn = 0.01 * u[3]
    ddn = 0.003 * u[4]
    dpd = 0.05 * u[2]
    dzd = 0.02 * u[3]

    return @SVector [dpn + dzn + ddn - dnp;
                     dnp - dpn - dpz - dpd;
                     dpz - dzn - dzd;
                     dpd + dzd - ddn]
end
u0_npzd = @SVector [8.0, 2.0, 1.0, 4.0]
"""
    prob_pds_npzd

Positive and conservative autonomous nonlinear PDS
```math
\\begin{aligned}
u_1' &= 0.01u_2 + 0.01u_3 + 0.003u_4 - \\frac{u_1u_2}{0.01 + u_1},\\\\
u_2' &= \\frac{u_1u_2}{0.01 + u_1}- 0.01u_2 - 0.5( 1 - e^{-1.21u_2^2})u_3 - 0.05u_2,\\\\
u_3' &= 0.5(1 - e^{-1.21u_2^2})u_3 - 0.01u_3 - 0.02u_3,\\\\
u_4' &= 0.05u_2 + 0.02u_3 - 0.003u_4
\\end{aligned}
```
with initial value ``\\mathbf{u}_0 = (8.0, 2.0, 1.0, 4.0)^T`` and time domain ``(0.0, 10.0)``.
There is one independent linear invariant, e.g. ``u_1+u_2+u_3+u_4 = 15.0``.

## References

- Hans Burchard, Eric Deleersnijder, and Andreas Meister.
  "Application of modified Patankar schemes to stiff biogeochemical models for the water column."
  Ocean Dynamics 55 (2005): 326-337.
  [DOI: 10.1007/s10236-005-0001-x](https://doi.org/10.1007/s10236-005-0001-x)
"""
prob_pds_npzd = ConservativePDSProblem(P_npzd, u0_npzd, (0.0, 10.0), std_rhs = f_npzd,
                                       linear_invariants = @SMatrix[1.0 1.0 1.0 1.0])

# stratospheric reaction problem
function P_stratreac(u, p, t)
    O1D, O, O3, O2, NO, NO2 = u

    Tr = 4.5
    Ts = 19.5
    T = mod(t / 3600, 24)
    if (Tr <= T) && (T <= Ts)
        Tfrac = (2 * T - Tr - Ts) / (Ts - Tr)
        sigma = 0.5 + 0.5 * cos(pi * abs(Tfrac) * Tfrac)
    else
        sigma = zero(t)
    end

    M = 8.120e16

    k1 = 2.643e-10 * sigma^3
    k2 = 8.018e-17
    k3 = 6.120e-4 * sigma
    k4 = 1.567e-15
    k5 = 1.070e-3 * sigma^2
    k6 = 7.110e-11
    k7 = 1.200e-10
    k8 = 6.062e-15
    k9 = 1.069e-11
    k10 = 1.289e-2 * sigma
    k11 = 1.0e-8

    r1 = k1 * O2
    r2 = k2 * O * O2
    r3 = k3 * O3
    r4 = k4 * O3 * O
    r5 = k5 * O3
    r6 = k6 * M * O1D
    r7 = k7 * O1D * O3
    r8 = k8 * O3 * NO
    r9 = k9 * NO2 * O
    r10 = k10 * NO2
    r11 = k11 * NO * O

    return @SMatrix [0.0 0.0 r5 0.0 0.0 0.0;
                     r6 r1+r10 r3 r1 0.0 0.0;
                     0.0 r2 0.0 0.0 0.0 0.0;
                     r7 r4+r9 r4+r7+r8 r3+r5 0.0 0.0;
                     0.0 0.0 0.0 0.0 0.0 r9+r10;
                     0.0 0.0 0.0 0.0 r8+r11 0.0]
end
function d_stratreac(u, p, t)
    O1D, O, O3, O2, NO, NO2 = u

    k2 = 8.018e-17
    k11 = 1.0e-8

    r2 = k2 * O * O2
    r11 = k11 * NO * O

    return @SVector [0.0, r11, 0.0, r2, 0.0, 0.0]
end
function f_stratreac(u, p, t)
    O1D, O, O3, O2, NO, NO2 = u

    Tr = 4.5
    Ts = 19.5
    T = mod(t / 3600, 24)
    if (Tr <= T) && (T <= Ts)
        Tfrac = (2 * T - Tr - Ts) / (Ts - Tr)
        sigma = 0.5 + 0.5 * cos(pi * abs(Tfrac) * Tfrac)
    else
        sigma = zero(t)
    end

    M = 8.120e16

    k1 = 2.643e-10 * sigma^3
    k2 = 8.018e-17
    k3 = 6.120e-4 * sigma
    k4 = 1.567e-15
    k5 = 1.070e-3 * sigma^2
    k6 = 7.110e-11
    k7 = 1.200e-10
    k8 = 6.062e-15
    k9 = 1.069e-11
    k10 = 1.289e-2 * sigma
    k11 = 1.0e-8

    r1 = k1 * O2
    r2 = k2 * O * O2
    r3 = k3 * O3
    r4 = k4 * O3 * O
    r5 = k5 * O3
    r6 = k6 * M * O1D
    r7 = k7 * O1D * O3
    r8 = k8 * O3 * NO
    r9 = k9 * NO2 * O
    r10 = k10 * NO2
    r11 = k11 * NO * O

    return @SVector [r5 - r6 - r7;
                     2 * r1 - r2 + r3 - r4 + r6 - r9 + r10 - r11;
                     r2 - r3 - r4 - r5 - r7 - r8;
                     -r1 - r2 + r3 + 2 * r4 + r5 + 2 * r7 + r8 + r9;
                     -r8 + r9 + r10 - r11;
                     r8 - r9 - r10 + r11]
end
u0_stratreac = @SVector [9.906e1, 6.624e8, 5.326e11, 1.697e16, 4e6, 1.093e9]
"""
    prob_pds_stratreac

Positive and nonconservative autonomous nonlinear PDS
```math
\\begin{aligned}
u_1' &= r_5 - r_6 -  r_7,\\\\
u_2' &= 2r_1 - r_2 + r_3 - r_4 + r_6 - r_9 + r_{10} - r_{11},\\\\
u_3' &= r_2 - r_3 - r_4 - r_5 - r_7 - r_8,\\\\
u_4' &= -r_1 -r_2 + r_3 + 2r_4+r_5+2r_7+r_8+r_9,\\\\
u_5' &= -r_8+r_9+r_{10}-r_{11},\\\\
u_6' &= r_8-r_9-r_{10}+r_{11},
\\end{aligned}
```
with reaction rates
```math
\\begin{aligned}
r_1 &=2.643⋅ 10^{-10}σ^3 u_4, & r_2 &=8.018⋅10^{-17}u_2 u_4 , & r_3 &=6.12⋅10^{-4}σ u_3,\\\\
r_4 &=1.567⋅10^{-15}u_3 u_2 , & r_5 &= 1.07⋅ 10^{-3}σ^2u_3,  & r_6 &= 7.11⋅10^{-11}⋅ 8.12⋅10^6 u_1,\\\\
r_7 &= 1.2⋅10^{-10}u_1 u_3, & r_8 &= 6.062⋅10^{-15}u_3 u_5, & r_9 &= 1.069⋅10^{-11}u_6 u_2,\\\\
r_{10} &= 1.289⋅10^{-2}σ u_6, & r_{11} &= 10^{-8}u_5 u_2,
\\end{aligned}
```
where
```math
\\begin{aligned}
T &= t/3600 \\mod 24,\\quad T_r=4.5,\\quad T_s = 19.5,\\\\
σ(T) &= \\begin{cases}1, & T_r≤ T≤ T_s,\\\\0, & \\text{otherwise}.\\end{cases}
\\end{aligned}
```

The initial value is ``\\mathbf{u}_0 = (9.906⋅10^1, 6.624⋅10^8, 5.326⋅10^{11}, 1.697⋅10^{16}, 4⋅10^6, 1.093⋅10^9)^T`` and the time domain ``(4.32⋅ 10^{4}, 3.024⋅10^5)``.
There are two independent linear invariants, e.g. ``u_1+u_2+3u_3+2u_4+u_5+2u_6=(1,1,3,2,1,2)\\cdot\\mathbf{u}_0`` and ``u_5+u_6 = 1.097⋅10^9``.

## References

- Stephan Nüsslein, Hendrik Ranocha, and David I. Ketcheson.
  "Positivity-preserving adaptive Runge-Kutta methods."
  Communications in Applied Mathematics and Computer Science 16 (2021): 155-179.
  [DOI: 10.2140/camcos.2021.16.155](https://doi.org/10.2140/camcos.2021.16.155)
"""
prob_pds_stratreac = PDSProblem(P_stratreac, d_stratreac, u0_stratreac, (4.32e4, 3.024e5),
                                std_rhs = f_stratreac,
                                linear_invariants = @SMatrix[1.0 1.0 3.0 2.0 1.0 2.0;
                                                             0.0 0.0 0.0 0.0 1.0 1.0])

function f_stratreac_scaled(u, p, t)
    uc = [9.906e1, 6.624e8, 5.326e11, 1.697e16, 4e6, 1.093e9]

    Tr = 4.5
    Ts = 19.5
    T = mod(t / 3600, 24)
    if (Tr <= T) && (T <= Ts)
        Tfrac = (2 * T - Tr - Ts) / (Ts - Tr)
        sigma = 0.5 + 0.5 * cos(pi * abs(Tfrac) * Tfrac)
    else
        sigma = zero(t)
    end

    M = 8.120e16

    k1 = 2.643e-10 * sigma^3
    k2 = 8.018e-17
    k3 = 6.120e-4 * sigma
    k4 = 1.567e-15
    k5 = 1.070e-3 * sigma^2
    k6 = 7.110e-11
    k7 = 1.200e-10
    k8 = 6.062e-15
    k9 = 1.069e-11
    k10 = 1.289e-2 * sigma
    k11 = 1.0e-8

    r1 = k1 * u[4] * uc[4]
    r2 = k2 * u[2] * uc[2] * u[4] * uc[4]
    r3 = k3 * u[3] * uc[3]
    r4 = k4 * u[3] * uc[3] * u[2] * uc[2]
    r5 = k5 * u[3] * uc[3]
    r6 = k6 * M * u[1] * uc[1]
    r7 = k7 * u[1] * uc[1] * u[3] * uc[3]
    r8 = k8 * u[3] * uc[3] * u[5] * uc[5]
    r9 = k9 * u[6] * uc[6] * u[2] * uc[2]
    r10 = k10 * u[6] * uc[6]
    r11 = k11 * u[5] * uc[5] * u[2] * uc[2]

    return @SVector [(r5 - r6 - r7) / uc[1];
                     (2 * r1 - r2 + r3 - r4 + r6 - r9 + r10 - r11) / uc[2];
                     (r2 - r3 - r4 - r5 - r7 - r8) / uc[3];
                     (-r1 - r2 + r3 + 2 * r4 + r5 + 2 * r7 + r8 + r9) / uc[4];
                     (-r8 + r9 + r10 - r11) / uc[5];
                     (r8 - r9 - r10 + r11) / uc[6]]
end
u0_stratreac_scaled = @SVector ones(6)
"""
    prob_ode_stratreac_scaled

Scaled version of the stratosperic reaction problem [`prob_pds_stratreac`](@ref).
Each component is scaled by its corresponding original initial value.  

The initial value is ``\\mathbf{u}_0 = (1,1,1,1,1,1)^T`` and the time domain ``(4.32⋅10^{4}, 3.024⋅10^5)``.

There are two independent linear invariants. The function `linear_invariants_stratreac_scaled` returns the invariance matrix.
"""
prob_ode_stratreac_scaled = ODEProblem(f_stratreac_scaled, u0_stratreac_scaled,
                                       (4.32e4, 3.024e5))

function linear_invariants_stratreac_scaled()
    return @SMatrix [99.06 6.624e8 1.5978e12 3.394e16 4.0e6 2.186e9;
                     0.0 0.0 0.0 0.0 4.0e6 1.093e9]
end

# mapk problem
function f_minmapk(u, p, t)
    k1 = 100 / 3
    k2 = 1 / 3
    k3 = 50
    k4 = 0.5
    k5 = 10 / 3
    k6 = 0.1
    k7 = 0.7

    return @SVector [k6 * u[6] - k7 * u[1] - k1 * u[1] * u[2] + k2 * u[4];
                     k5 * u[3] - k1 * u[1] * u[2];
                     k2 * u[4] - k3 * u[1] * u[3] + k4 * u[5] - k5 * u[3];
                     k1 * u[1] * u[2] - k2 * u[4];
                     k3 * u[1] * u[3] - k4 * u[5];
                     k7 * u[1] - k6 * u[6]]
end
function P_minmapk(u, p, t)
    k1 = 100 / 3
    k2 = 1 / 3
    k3 = 50
    k4 = 0.5
    k5 = 10 / 3
    k6 = 0.1
    k7 = 0.7

    return @SMatrix [0 0 0 k2*u[4] 0 k6*u[6]
                     0 0 k5*u[3] 0 0 0;
                     0 0 k2*u[4] 0 k4*u[5] 0;
                     k1*u[1]*u[2] 0 0 0 0 0;
                     0 0 k3*u[1]*u[3] 0 0 0;
                     k7*u[1] 0 0 0 0 0]
end
function D_minmapk(u, p, t)
    k1 = 100 / 3

    return @SVector [0; k1 * u[1] * u[2]; 0; 0; 0; 0]
end

u0 = @SVector [0.1; 0.175; 0.15; 1.15; 0.81; 0.5]
tspan = (0.0, 200.0)

"""
    prob_pds_minmapk

Positive and nonconservative autonomous nonlinear PDS
```math
\\begin{aligned}
u_1' &= k_6u_6-k_7u_1-k_1u_1u_2 +k_2u_4,\\\\
u_2' &= k_5u_3-k_1u_2u_2,\\\\
u_3' &= k_2u_4-k_3u_1u_3+k_4u_5-k_5u_3,\\\\
u_4' &= k_1u_1u_2-k_2u_4,\\\\
u_5' &= k_3u_1u_3-k_4u_5,\\\\
u_6' &= k_7u_1-k_6u_6,
\\end{aligned}
```
with constants
```math
\\begin{aligned}
k_1 &=\\frac{100}{3}, & k_2 &=\\frac{1}{3}, & k_3 &=50,\\\\
k_4 &=0.5, & k_5 &=\\frac{10}{3} ,  & k_6 &= 0.1,\\\\
k_7 &= 0.1.
\\end{aligned}
```

The initial value is ``\\mathbf{u}_0 = (0.1, 0.175, 0.15, 1.15, 0.81, 0.5)^T`` and the time domain ``(0, 200)``.
There are two independent linear invariants, e.g. ``u_1+u_4+u_6=1.75`` and ``u_2+u_3+u_4+u_5 =2.285``.

## References

- Sergio Blanes, Arieh Iserles, and Shev Macnamara.
  "Positivity preserving methods for ordinary differential equations."
  ESAIM: Mathematical Modelling and Numerical Analysis 56 (2022): 1843–1870.
  [DOI: 10.1051/m2an/2022042](https://doi.org/10.1051/m2an/2022042)
- Otto Hadač, František Muzika, Vladislav Nevoral, Michal Přibyl, and Igor Schreiber
  "Minimal oscillating subnetwork in the Huang-Ferrell model of the MAPK cascade."
  PLoS ONE 12 (2017): e0178457.
  [DOI: 10.1371/journal.pone.0178457](https://doi.org/10.1371/journal.pone.0178457)
"""
prob_pds_minmapk = PDSProblem(P_minmapk, D_minmapk, u0, tspan; std_rhs = f_minmapk,
                              linear_invariants = @SMatrix[1.0 0.0 0.0 1.0 0.0 1.0;
                                                           0.0 1.0 1.0 1.0 1.0 0.0])

# SIRD (Sen & Sen) problem

function f_sird_sensen(u, p, t)
    Npop = 6.046e7
    alpha = 0.0194
    beta = 7.567
    mu = 2.278e-6
    eta = 9.180e-7
    sigma = 1.4633e-3
    tau = 1.109e-4
    xi = 0.263
    gamma = 0.021
    delta = 0.077
    lambda = 6.2800e-04
    Kd = 0.0013

    # infection-like term
    inf_term = (beta * u[5] + sigma * u[2]) / Npop + eta

    return @SVector [-alpha * u[1] - inf_term * u[1];
                     xi * u[4] - tau * u[2];
                     alpha * u[1] - mu * u[3];
                     inf_term * u[1] + mu * u[3] - (gamma + xi) * u[4];
                     tau * u[2] + gamma * u[4] - delta * u[5];
                     lambda * u[7];
                     delta * u[5] - (lambda + Kd) * u[7];
                     Kd * u[7]]
end

function P_sird_sensen(u, p, t)
    Npop = 6.046e7
    alpha = 0.0194
    beta = 7.567
    mu = 2.278e-6
    eta = 9.180e-7
    sigma = 1.4633e-3
    tau = 1.109e-4
    xi = 0.263
    gamma = 0.021
    delta = 0.077
    lambda = 6.2800e-04
    Kd = 0.0013

    # P[i,j] is flux from compartment j -> i
    P41 = (u[1] / Npop) * (beta * u[5] + sigma * u[2]) + eta * u[1]

    return @SMatrix [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
                     0.0 0.0 0.0 xi*u[4] 0.0 0.0 0.0 0.0;
                     alpha*u[1] 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
                     P41 0.0 mu*u[3] 0.0 0.0 0.0 0.0 0.0;
                     0.0 tau*u[2] 0.0 gamma*u[4] 0.0 0.0 0.0 0.0;
                     0.0 0.0 0.0 0.0 0.0 0.0 lambda*u[7] 0.0;
                     0.0 0.0 0.0 0.0 delta*u[5] 0.0 0.0 0.0;
                     0.0 0.0 0.0 0.0 0.0 0.0 Kd*u[7] 0.0]
end

# initial value (from SirdTest with sost = 1e-10)
u0_sird_sensen = @SVector [6.046e7 - (4e-10 + 3.0); 1e-10; 1e-10; 1.0; 1.0; 1e-10; 1.0;
                           1e-10]

tspan_sird_sensen = (0.0, 180.0)

"""
    prob_pds_sird_sensen

Positive and conservative autonomous nonlinear PDS
```math
\\begin{aligned}
u_1' &= -\\alpha u_1 - \\left(\\frac{\\beta}{N_{\\mathrm{pop}}}u_5 + \\frac{\\sigma}{N_{\\mathrm{pop}}}u_2 + \\eta\\right)u_1,\\\\
u_2' &= \\xi u_4 - \\tau u_2,\\\\
u_3' &= \\alpha u_1 - \\mu u_3,\\\\
u_4' &= \\left(\\frac{\\beta}{N_{\\mathrm{pop}}}u_5 + \\frac{\\sigma}{N_{\\mathrm{pop}}}u_2 + \\eta\\right)u_1
        + \\mu u_3 - (\\gamma + \\xi)u_4,\\\\
u_5' &= \\tau u_2 + \\gamma u_4 - \\delta u_5,\\\\
u_6' &= \\lambda u_7,\\\\
u_7' &= \\delta u_5 - (\\lambda + K_d)u_7,\\\\
u_8' &= K_d u_7,
\\end{aligned}
```

with constants

```math
\\begin{aligned}
N_{\\mathrm{pop}} &= 6.046\\cdot 10^{7}, \\\\
\\alpha &= 0.0194, \\\\
\\beta  &= 7.567, \\\\
\\mu    &= 2.278\\cdot 10^{-6}, \\\\
\\eta   &= 9.180\\cdot 10^{-7}, \\\\
\\sigma &= 1.4633\\cdot 10^{-3}, \\\\
\\tau   &= 1.109\\cdot 10^{-4}, \\\\
\\xi    &= 0.263, \\\\
\\gamma &= 0.021, \\\\
\\delta &= 0.077, \\\\
\\lambda &= 6.2800\\cdot 10^{-4}, \\\\
K_d    &= 0.0013.
\\end{aligned}
```

The initial value is ``\\mathbf{u}_0 = (6.046\\cdot 10^{7}-(4\\cdot 10^{-10}+3),\\,10^{-10},\\,10^{-10},\\,1,\\,1,\\,10^{-10},\\,1,\\,10^{-10})^T`` and the time domain ``(0.0, 180.0)``.
There is one independent linear invariant, namely total population ``u_1+u_2+u_3+u_4+u_5+u_6+u_7+u_8 = N_{\\mathrm{pop}}``.

## References

- D. Sen and D. Sen.
  "Use of a modified SIRD model to analyze COVID-19 data."
  *Industrial & Engineering Chemistry Research* 60(11) (2021): 4251–4260.
  [DOI: 10.1021/acs.iecr.0c04754](https://doi.org/10.1021/acs.iecr.0c04754) :contentReference[oaicite:0]{index=0}

- Giuseppe Izzo, Eleonora Messina, Mario Pezzella, and Antonia Vecchio.
  "Modified Patankar Linear Multistep Methods for Production-Destruction Systems."
  *Journal of Scientific Computing* 102 (2025): 87.
  [DOI: 10.1007/s10915-025-02804-5](https://doi.org/10.1007/s10915-025-02804-5) :contentReference[oaicite:1]{index=1}
"""

prob_pds_sird_sensen = ConservativePDSProblem(P_sird_sensen, u0_sird_sensen,
                                              tspan_sird_sensen, std_rhs = f_sird_sensen,
                                              linear_invariants = @SMatrix[1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0])

# diffusion problem
function f_diffusion!(du, u, p, t)
    dx = p.dx
    K = p.K_ev
    N = length(u)
    invdx2 = one(eltype(u)) / (dx^2)

    @inbounds begin
        du[1] = (K[2] * u[2] - K[1] * u[1]) * invdx2
        for i in 2:(N - 1)
            du[i] = (K[i + 1] * u[i + 1] +
                     K[i - 1] * u[i - 1] -
                     2 * K[i] * u[i]) * invdx2
        end
        du[N] = (K[N - 1] * u[N - 1] - K[N] * u[N]) * invdx2
    end
    return nothing
end

function P_diffusion!(P::Tridiagonal, u, p, t)
    dx = p.dx
    K = p.K_ev
    N = length(u)
    invdx2 = one(eltype(u)) / (dx^2)

    fill!(P.dl, zero(eltype(P)))
    fill!(P.d, zero(eltype(P)))
    fill!(P.du, zero(eltype(P)))

    @inbounds for i in 1:(N - 1)
        P.du[i] = K[i + 1] * u[i + 1] * invdx2
        P.dl[i] = K[i] * u[i] * invdx2
    end

    return nothing
end

N_diffusion = 100
L_diffusion = 1.0
dx_diffusion = L_diffusion / N_diffusion
x_diffusion = dx_diffusion / 2 .* ones(N_diffusion)
for j in 2:N_diffusion
    x_diffusion[j] = x_diffusion[j - 1] + dx_diffusion
end

D0 = 1e-2
kfun = x -> 1e-5 +
            (x - 2 * L_diffusion / 3)^2 * D0 *
            atan(0.5 * (2 * x - 3 * L_diffusion)) /
            (0.5 * (2 * x - 3 * L_diffusion))
K_ev_diffusion = kfun.(x_diffusion)

f0 = x -> 2 * (1 - sin(pi * (x * pi / 2 - 0.25))^2)
u0_diffusion = [f0(xi) for xi in x_diffusion]

tspan_diffusion = (0.0, 3.0)

p_diffusion = (dx = dx_diffusion, K_ev = K_ev_diffusion)

p_prototype_diffusion = Tridiagonal(zeros(eltype(u0_diffusion), N_diffusion - 1),
                                    zeros(eltype(u0_diffusion), N_diffusion),
                                    zeros(eltype(u0_diffusion), N_diffusion - 1))

"""
    prob_pds_diffusion

Positive and conservative autonomous nonlinear production–destruction system
obtained from a finite-volume discretization of a one-dimensional diffusion equation
with spatially varying diffusion coefficient.

```math
\\begin{aligned}
u_i' &= \\sum_{j=1}^{N} \\bigl( P_{ij}(u) - P_{ji}(u) \\bigr), \\qquad i = 1,\\dots,N,\\\\
P_{i,i+1}(u) &= \\frac{1}{\\Delta x^2} K_{i+1} u_{i+1},\\qquad
P_{i+1,i}(u) = \\frac{1}{\\Delta x^2} K_i u_i,
\\end{aligned}
```

with ``P_{i,j}(u)=0`` otherwise.


The grid consists of N = 100 cells with width ``\\Delta x = 10^{-2}``
and centers ``x_i = (i-\\tfrac12)\\Delta x`` (``L = 1``).
The initial value is ``\\mathbf{u}_0 = (u_1^0,\\dots,u_N^0)^T`` with
``u_i^0 = f(x_i)``, and the time domain ``(0.0, 3.0)``.

There is one independent linear invariant, namely
``\\sum_{i=1}^{N} u_i = \\text{const}.``

## References

- Giuseppe Izzo, Eleonora Messina, Mario Pezzella, and Antonia Vecchio.
  "Modified Patankar Linear Multistep Methods for Production-Destruction Systems."
  *Journal of Scientific Computing* 102 (2025): 87.
  [DOI: 10.1007/s10915-025-02804-5](https://doi.org/10.1007/s10915-025-02804-5) :contentReference[oaicite:1]{index=1}
"""
prob_pds_diffusion = ConservativePDSProblem(P_diffusion!,
                                            u0_diffusion,
                                            tspan_diffusion,
                                            p_diffusion;
                                            p_prototype = p_prototype_diffusion,
                                            std_rhs = f_diffusion!,
                                            linear_invariants = ones(1, N_diffusion))