# Spectral Collocation (or Weighted Residual) Methods to solve Stiff ODEs?

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.

• Could you specify precisely what you mean by "projection methods", with equations or a reference? – David Ketcheson Nov 4 '15 at 7:15
• @DavidKetcheson See my edits. Thanks for looking at this! – jlperla Nov 4 '15 at 16:58
• I can't tell from what you have written whether you have an initial or boundary value problem. Since you refer to stiffness and Runge-Kutta methods, I would guess it is an initial value problem. But if so, this is a very expensive method since it couples all the time steps into one big system. – David Ketcheson Nov 4 '15 at 19:32
• By the way, "collocation" and "projection" refer to very many different things in the realm of numerical methods. In my view, what you are proposing is quite different from what is explained on the Wikipedia page you have linked to. – David Ketcheson Nov 4 '15 at 19:35
• The wikipedia article is saying exactly what I mean: "choose a finite-dimensional space of candidate solutions (usually, polynomials up to a certain degree) and a number of points in the domain (called collocation points), and to select that solution which satisfies the given equation at the collocation points." – jlperla Nov 4 '15 at 22:23