# Solving an ODE while maintaining weak positivity and weak monotonicity

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.

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.

• Thank you so much. Are there any tricks or best practices when using spectral methods? Nov 8 '15 at 15:14
• I don't know of any spectral methods where it is possible to ensure positivity for any general class of problems. Nov 9 '15 at 16:28
• Thank you so much. I am fine with that accuracy. For finite differences: the final strange part of this problem as listed is the asymptotic Convergence above, where $F(\infty)$ is a constant. It is great to hear that backwards euler preserves this, but part of the problem I have run into before is that since $\lim{z\to\infty}F'(z) = 0$, you have what I believe is called a singular point(?) in the ODE. this has caused havoc for large $z$ as the ODE becomes asymptotically flat. Is Backwards Euler stable for this issue (singular point exactly at the boundary we are working back from?) Nov 9 '15 at 17:50
• Now you say you are working back from the boundary, but in the question you seem to be solving forward from zero. Nov 9 '15 at 18:33
• Sorry, I was responding to your comment that `Indeed, if the exact solution satisfies F(z)≤1, Backward Euler will preserve that too. ' My comment was what if that boundary is at infinity (choosing a large max z), so that F'(z) -> 0 for a large segment of whatever domain I give it. This has cause significant problems with other methods. (FYI: as I know the analytical boundaries, I could set it up as a BVP, but I am not sure that helps anything here) Nov 9 '15 at 18:38