I have a system of ODEs which is (at least moderately) stiff.
Consider the class of spectral collocation methods https://en.wikipedia.org/wiki/Spectral_method or the related class of weighted residual methods.
Basically, take a system of ODEs like $$ y'(t) = f(t, y(t)) $$ where $f(\cdot,\cdot)$ is nonlinear in my case and leads to a stiff ODE, and I want to solve it on $t \in [0,M]$. Then approximate the solution $$y(t) \approx \sum_{i=0}^N a_i T_i(t)$$ where $T_i(\cdot)$ are a basis (e.g. Chebyshev polynomials), $N$ is the maximum order, and the basis is over $[0,M]$ if appropriate. Note that given $\{a_i\}$, $y'(t)$ can be calculated by differentiating the polynomial basis.
Now, evaluate the ODE given a $\{a_i\}$ at the roots of the polynomial basis (which can be shown to be optimal) and find the residual. Choose $\{a_i\}$ to minimize this residual, or just solve a system of equations with the residual = 0.
Question: how well does spectral collocation with a Chebyshev basis handle systems that seem stiff when solved with finite differences? The worry is that the same issues that cause finite differences to have issues (e.g. convergence to a very flat function with no slope) could cause major issues with a finite dimensional approximation of the function.
If Chebyshev has problems, then is there a better basis which might have less issues.
A little more context: I have implemented spectral collocation with Chebyshev polynomials to solve a stationary Kolmogorov Forward equation in CDFs. These, of course, need to be strictly positive and weakly monotonically increasing and asymptotically go towards $1$ if it is a proper probability distribution. I am seeing weird behavior where the pdf drops below $0$ (i.e., it becomes decreasing) or strange jumps at the corners. It is very hard for me to tell if these are artifacts in the solution scheme or if I have just chosen parameters for the model where no solution exists. Increasing the number of Chebyshev basis functions doesn't eliminate them, but that could be because they can't, or only slowly, converge.
Even more context: For what it is worth, in the simplest form of my problem there are two sets of ODEs (which are ultimately coupled):
Parameters: $\lambda_{\ell}, \lambda_h, r,$ and $\chi$. In some ways $F_{\ell}'(0), F_h'(0),$ and $g$ are also parameters for the purposes here. There are two discrete states which leads to the system: $i \in \{\ell,h\}$
The variable $z\in[0,\infty)$ in reality, but these converge pretty quickly so choosing some $z \in [0,\bar{z}]$ for collocation methods is reasonable.
System (1) The following comes from a stationary Kolmogorov Forward Equation (given a $\gamma(z)$ function such that $\gamma(0) = 0$ and $\lim\limits_{z\to\infty} \gamma(z) = g$. The latter condition leads to a vanishing derivative term in the ODE (i.e., a singular mass matrix if put in that canonical form.) $$ \begin{align} 0 &= g F_{\ell}'(z) + \lambda_h F_h(z) - \lambda_{\ell}F_{\ell}(z) + g (F_{\ell}'(0) + F_h'(0))(F_{\ell}(z) + F_h(z)) - g F_{\ell}'(0)\\ 0 &= (g - \gamma(z))F_h'(z) + \lambda_{\ell}F_{\ell}(z) - \lambda_h F_h(z) - g F_h'(0) \end{align} $$
- Forget that this is a linear ODE, as the actual one is nonlinear and much trickier. However, it can still be written in a $M(z, F(z))\cdot F'(z) = \Phi(z, F(z))$ term for some operators and a mass matrix.
- The trouble seems to come out of the term on the derivative going to $0$ (i.e., a singular mass matrix).
- It might seem odd that there is a $F_{\ell}'(0)$ and a $F_h'(0)$ in the ODE as parameters but this is not a mistake. You can easily show that any solution to the ODE will have the $F_i'(0)$ matching these constants by construction.
- The initial condition is $F_i(0) = 0$. Furthermore, from any $F_i'(0)$, one can analytically find $\lim\limits_{z\to\infty}F_h(z)$. With $F_{\ell}(\infty) + F_h(\infty) = 1$, this means we can write it as a BVP if we wish.
System (2) The following comes from a Hamilton-Jacobi-Bellman equation
\begin{align} 0 &= 1 - (r + \lambda_{\ell})w_{\ell}(z) - g w'_{\ell}(z) + \lambda_h w_h(z)\\ 0 &= 1 - (r + \lambda_h)w_h(z) - (g - \frac{\chi}{2} w_h(z))w_h'(z) + \lambda_h w_{\ell}(z) + \frac{\chi}{4} w_h(z)^2 \end{align} subject to, $w_{\ell}(0) = w_h(0) = 0$ and defined on $z \in [0,\infty)$. You can analytically find $w_{\ell}(\infty)$ and $w_h(\infty)$ which can help with various methods to set it as a BVP.
One can show that $\lim\limits_{z\to\infty}(g - \frac{\chi}{2} w_h(z)) = 0$ which leads to a vanishing derivative term in the ODE.
Coupled: In reality, the $\gamma(z)$ from the first solution is $$ \gamma(z) \equiv \frac{\chi}{2}w_h(z) $$
Leaving Out: A whole bunch of equilibrium conditions in terms of integrals of $w_i(z)$ and $F_i(z)$ which pin down $g, F'_{\ell}(0)$ and $F'_h(0)$. Evaluating these concurrently with the solution to the ODEs is the reason that spectral methods are preferable since I can add in more equations and solve everything at the same time.
... is what I have written useful? Maybe not, but you can see the asymptotic singularity in the mass matrix. Is this stiff? I think so. When you solve the second system with finite differences, for example, it has all the hallmarks of a stiff system.
