Consider an ODE of the form

$$c_0 u + c_1 u' + c_2 u'' \enspace = \enspace f(x,u) \quad .$$

I want to solve this ODE with spectral methods. To this end, let $\{ x_k \}_{k=0}^N$ be suitable nodes, e.g. Tchebychev-Radau points $x_k = \cos\big( \tfrac{2\pi k}{2N+1} \big)$. I approximate $u$ by a series of the form

$$ u(x) \enspace = \enspace \sum_{i=0}^N u_i L_i(x) $$

with $u_i = u(x_i)$ and where $L_i(x)$ are Lagrange-Polynomials with nodes $x_k$ such that $L_i(x_k) = \delta_{ik}$.

1.)$\enspace$ Consider the linear case $f(x,u) \equiv 0$. Evaluation at the nodes $x_k$ allows rewriting the ODE as a linear system of equations

$$\sum_{j=0}^N \Big( c_0 \delta_{ij} + c_1 D_{ij} + c_2 D^2_{ij} \Big) u_j \enspace = \enspace 0 \quad , \qquad i \in \{0,\ldots, N\}$$

with $D_{ij}$ being the differentiation matrix. With given inital conditions, this can be solved.

2.)$\enspace$ However, assume now some nonlinearity, e.g. $f(x,u) = u^3$, then this ODE can be written as

$$ \sum_{j=0}^N \Big[ c_0 \delta_{ij} + c_1 D_{ij} + c_2 D_{ij}^2 \Big] u_j \enspace = \enspace u_i^3 $$

How can this ODE be solved with (pseudo-)spectral (collocation) methods? Due to the nonlinearity, inverting the matrix on the LHS is no longer an option. Any ideas?


The rewritten ODE is obtained as follows: Insertion of the expansion into the ODE yields

$$ \sum_{k=0}^N \Big[ c_0 u_k + c_1 u_k' + c_2 u_k'' \Big] L_k(x) \enspace = \enspace \sum_{k=0}^N u_k^3 L_k(x)$$

where $u_i' = u'(x_i)$ and $u_i'' = u''(x_i)$. Evaluation at Tchebychev-Radau points $x_i$:

$$ c_0 u_i + c_1 u_i' + c_2 u_i'' \enspace = \enspace u_i^3 $$

since $L_k(x_i) = \delta_{ki}$. With the differentiation matrices

$$ \sum_{j=0}^N \Big[ c_0 \delta_{ij} + c_1 D_{ij} + c_2 D_{ij}^2 \Big] \enspace = \enspace u_i^3 $$

  • 1
    $\begingroup$ Unfortunately, I dont think the rewritten ODE is correct in the general case. I dont see how you could possible transform $f(x,\sum u_iL_i)\to\sum L_i f(x,u_i)$ for nonlinear $f$. Also, it is unclear what $x_i$ is. $\endgroup$
    – maxmilgram
    Sep 30 at 11:19
  • $\begingroup$ Maybe I don't get it, but your $u_i$ are not functions anymore, are they? That's the whole point of the approximation... $\endgroup$
    – maxmilgram
    Sep 30 at 18:01
  • $\begingroup$ I updated the question (again) in order to (hopefully) clear it up. $\endgroup$
    – Octavius
    Sep 30 at 18:17
  • $\begingroup$ Are there any particular boundary conditions you are interested in? Surely that will affect the well-posedness of any numerical discretization $\endgroup$
    – whpowell96
    Sep 30 at 19:18

2 Answers 2


Inserting the Ansatz into the differential equation yields: $$ \sum_{j=0}^N \Big[ c_0 \delta_{ij} + c_1 D_{ij} + c_2 D_{ij}^2 \Big] u_j \enspace = \enspace f(x_i,u_i) $$ Since your problem is nonlinear you also need a solver for nonlinear problems most notably and commonly used Newton's method:

Indeed by defining

$$ F_i(\mathbf{u}) = \sum_{j=0}^N \Big[ c_0 \delta_{ij} + c_1 D_{ij} + c_2 D_{ij}^2 \Big] u_j - f(x_i,u_i)\\ J_{ij}(\mathbf{u})= c_0 \delta_{ij} + c_1 D_{ij} + c_2 D_{ij}^2 - \delta_{ij}\frac{\partial}{\partial u}f(x_i,u_i) $$

We iteratively solve for $$ J(\mathbf{u}_{n})(\mathbf {u} _{n+1}-\mathbf {u} _{n})=-F(\mathbf {u} _{n}) $$

Depending on you inital conditions you will need to add these to $F$ (and $J$) e.g. for homogenous Neumann boundary conditions: $$ F_i(\mathbf{u})=\begin{cases}\sum_{i=0}^N u_i L'(0) =0\\ \sum_{i=0}^N u_i L'(1) =0\\ \sum_{j=0}^N \Big[ c_0 \delta_{ij} + c_1 D_{ij} + c_2 D_{ij}^2 \Big] u_j - f(x_i,u_i) \end{cases} $$ However the stability of the scheme will usually depend on the boundary conditions. I'd be happy to provide some numerical insights if you provide parameters $c_i$ and boundary conditions.


Choosing as the basis functions the Lagrange interpolants,

$$ l_i(x) = \prod_{j=0\\ j\ne i}^n\frac{x-x_j}{x_i-x_j} $$

we have

$$ u(x) \approx \sum_i u_i l_i(x) $$

so the ODE residual at point $x=x_i$ is

$$ \delta_k = c_0\sum_i u_i l_i(x_k)+c_1\sum_i u_i l_i'(x_k)+c_2\sum_i u_il_i''(x_k)-u_k^3 $$

Here $f(x,u) = u^3$ and $l_i^{(n)}(x_k)$ should be interpreted as first calculating the derivative $l_i^{(n)}(x)$ and then $x=x_k$. The boundary conditions $u(a) = u_a, u(b) = u_b$ can be handled as $u_0 = u_a, u_n = u_b$ and the other $u_k$ values can be obtained thru a minimization procedure

$$ u^* = \arg\min_{u}\left((u_0-u_a)^2+(u_n-u_b)^2+\sum_k\delta_k^2\right) $$


$$ \cases{ l_i'(x) = l_i(x)\sum_{k=0\\ k \ne i}^{k=n}\frac{1}{x-x_k}\\ l_i''(x) = l_i(x)\left(\left(\sum_{k=0\\ k \ne i}^{k=n}\frac{1}{x-x_k}\right)^2-\sum_{k=0\\ k \ne i}^{k=n}\frac{1}{(x-x_k)^2}\right) } $$


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .