Let us considered the following system of ODEs \begin{align*} \dfrac{dX}{dt} = f(X), \tag{1.1} \end{align*} where the unknown $X\in \mathcal{D} \subset \mathbb{R}^l$ and it is stiff. However, for simplicity, I would like to suppose that $X$ is a scalar function, i.e, $X(t) \in \mathbb{R}$ and $f$ has the same dimension as $X$. There is quite effective traditional way to solve this equation is the "Implicit Euler method", which solves the equation with the following scheme \begin{align*} &\dfrac{X^{n+1}-X^n}{\Delta t} = f\left(X^{n+1}\right)\\ \Leftrightarrow & X^{n+1} - X^n - \Delta t f(X^{n+1}) = 0 \end{align*} This method sets us free from any stability condition, which enables us to examine larger time steps in order to save computational costs. However, the finding of the result of the algebraic system in $X^{n+1}$ becomes harder when the large time steps are considered. The complexity is that we do not know whether or not a solution exists at all. However, the finding of the result of the algebraic system in $X^{n+1}$ becomes harder when the large time steps are considered. The complexity is that we do not know whether or not a solution exists at all. In the situation of existence, the Newton method is the most popular method for solving the discretized system, however, it may not converge or converge extremely slowly. In contrast, when there is no solution, we have to restart the current iteration by dividing the time step by two, for example. For ensuring the convergence of the Newton method, we need to choose a good time step, but this work is not easy. Therefore, I would like to propose another scheme based on the complementarity condition as follows \begin{align*} X^{n+1} - X^n - (1-v)\Delta t f(X^{n+1}) & = 0 \tag{1.2a}\\ 1 - (1-v)\Delta t f'(X^{n+1}) -\varepsilon - w &= 0 \in \mathbb{R} \tag{1.2b} \\ \min(v,w) &= 0 \in \mathbb{R} \tag{1.2c} \end{align*} where $v$ and $w$ are two new unknowns to be sought for simultaneously to $X^{n+1}$. The variables $v$ and $w$ are introduced to rewrite the system in an equivalent way that allows us to work with a possibly smaller time-step $(1 — v )\Delta t$ instead of $\Delta t$. The reason why we wish to reduce the time-step in this fashion is to guarantee that the actual time-step $(1 — v )\Delta t$ is the largest possible one for which the Jacobian matrix of the system remains invertible. To understand the last point, let us imagine that $v$ is fixed but $\Delta t$ is very small, or that is $\Delta t$ fixed but $v$ is very close to 1. Then, the derivative of (1.2a) with respect to X, which is the left-hand side of (1.2b), equal to $1 – (1 — v )\Delta t f'(X)$, is very close to 1. When $(1 — v )\Delta t$ moves away from 0, its value may vanish and become negative. We want to stop just before its value becomes 0, to avoid crossing the frontier where the system (1.4a) is singular. The idea of the complementarity condition (1.2c) is that we have basically two cases:
- If $v = 0$ and $w > 0$, then there is no need to reduce the time-step since the derivative of (1.2a) is already invertible.
- If $w = 0$ and $v > 0$, then it is necessary to reduce the time-step and the new system has reached the limit zone. Here, maybe a small safety coefficient $\varepsilon > 0$ could be useful and we’ll probably have to work with $\min(v, w—\varepsilon ) = 0$ for (1.2c).
I have tested this scheme with two existing methods are "The Newton min method" and the "Single-stage interior point method". The result seems to be good for the case nonlinear example $f(X) = -aX(1-X)$, where $a>0$. However, when I use this scheme for the case $f(X) = -aX$, it cannot deal with the stiffness as well as I imagine. My code is in
https://drive.google.com/file/d/19N7puarSUVS-KDPJ2k6_3q2485cDK8EZ/view?usp=share_link
and the model I have used in this code is model = 3. I hope that there are some opinions and ideas to improve this scheme a little bit. In addition, there is still no theoretical proof of the efficiency of this scheme for the moment.
