# Nonlinear ODE with (pseudo-)spectral method

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

• 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. Sep 30 at 11:19
• Maybe I don't get it, but your $u_i$ are not functions anymore, are they? That's the whole point of the approximation... Sep 30 at 18:01
• I updated the question (again) in order to (hopefully) clear it up. Sep 30 at 18:17
• Are there any particular boundary conditions you are interested in? Surely that will affect the well-posedness of any numerical discretization Sep 30 at 19:18

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

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