Nonlinear ODE with (pseudo-)spectral method

Question

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?

---

DETAILS:

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 $$

Answer 1

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.

Answer 2

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) $$

NOTE

$$ \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) } $$