Solving an ODE while maintaining weak positivity and weak monotonicity

Question

I have a system of $N$ ODEs of the form, $$ M(z,F(z)) \cdot F'(z) = \Phi(z,F(z)) $$ where the mass matrix is $M(z,F): R\times R^N \to R^{N\times N}$ and $\Phi(z,F):R\times R^N \to R^N$ is (potentially) nonlinear. The initial condition is $F(0) = 0_N$. I need to solve this for an arbitrarily large $z$, but it is stationary so I can choose some $z_{\max}$.

A few notes:

• $M(z,F)$ is diagonal but it is asymptotically singular (i.e., as $z\to\infty$, a few of the diagonals approach $0$). This is why I left it in the mass matrix form, but I could just as easily write it as a $F'(z) = \Phi(z,F)$ form by dividing by the diagonal.
• This comes from a stationary Kolmogorov Forward equation in CDFs $F(z)$ (where the $N$ equations are an additional discrete state in the distribution, but I don't think that really matters here). This leads to a few crucial requirements on the solution:
  • $F(0) = 0_N$ since the minimum of support of the distribution is at $0$
  • Weak Positivity: $F(z) \geq 0$ for all $z$ to be a valid CDF, where the only $0$ in practice will be at the initial condition.
  • Weak Monotonicity: $F'(z) \geq 0$ for all $z$, since PDFs can't be negative.
  • Convergence: $\lim_{z\to\infty}F(z) < \infty$ to ensure that this is a valid CDF. (In reality, $\lim_{z\to\infty}\sum_{n=1}^{N}F_n(z) = 1$, but I can deal with this through other means)

Is this an IVP or a BVP? Without explaining the details of the setup, the short answer is that I can choose whichever is more convenient and robust.

QUESTION: What is the best numerical algorithm (or transformation of the ODE) to solve this problem? In particular, I would love to maintain positivity and monotonicity, and I suspect that the convergence requirement is making the problem stiff as $z$ gets large (which could also be related to the singular mass matrix).

I am willing to sacrifice speed (and even accuracy) for robustness here, but would prefer to use an algorithm with off-the-shelf software.

Answer 1

I assume that the exact solution satisfies the positivity and monotonicity conditions.

The most robust method that will guarantee these properties is the Backward Euler method. It will maintain positivity and monotonicity under arbitrarily large step sizes. Indeed, if the exact solution satisfies $F(z)\le 1$, Backward Euler will preserve that too. It is only first-order accurate, so you may need to take small step sizes anyway for the sake of accuracy.

For arbitrary nonlinear $\Phi$, it is known that there exists no 2nd-order method that preserves positivity under arbitrary step sizes.

If you can show that positivity is preserved under a forward Euler step (with some upper limit on the step size, obviously) then the class of Strong stability preserving (SSP) methods can be used to get higher-order accuracy, but you will have to deal with a restriction on the time step size, and you are likely to lose the high order accuracy anyway due to the asymptotic singularity.