What is the correct way to implement Neumann type boundary conditions for solving a PDE using Chebyshev's collocation method?

Question

I am studying spectral methods for solving PDEs numerically. I finished a chapter that explains how to use Chebyshev's collocation method to solve them. Though the explanation in the book is quite good, I still have some questions regarding the implementation of Neumann type boundary conditions (both homogeneous and non-homogeneous type). I would like to ask here for some help on clarifying the problem.

It would probably be best to use a specific example. Let's say we would like to solve the heat equation:

$$T_t = T_{xx} \tag 1$$

with a zero initial condition and the following boundary conditions:

$$T_x(-1,t) = f(t) \tag 2$$ $$T_x(1,t) = 0 \tag 3$$

Using the Euler's method to solve for time and Chebyshev's differentiation matrix $D$ for spatial derivatives, we simply write:

$$ T_{i+1}=T_i+\Delta tD^2T_i \tag 4$$

Based on what I read, I thought the boundary conditions can be incorporated the following way. I would use equation $(4)$ to calculate the temperature at the next time step and then modify the first and last element of the vector $T_{i+1}$ before proceeding to the next step. I thought I should do the modification by setting:

$$T_{i+1}(1) = f_i-D(1, :)T_{i+1}(1) \tag 5$$ $$T_{i+1}(\mathrm{end}) = 0-D(\mathrm{end}, :)T_{i+1}(\mathrm{end}) \tag 6$$

I used Matlab notation. The $D(1, :)$ means all columns of the first row and $D(\mathrm{end}, :) $ means all columns of the last row of the Chebyshev's differentiation matrix. However, when I noticed the dimensions were incorrect in $(5)$ and $(6)$ for matrix multiplication, I realized that this is incorrect. So, I don't understand how to implement these conditions.

What is the correct way to implement these boundary conditions into the solution algorithm, and why? A short derivation with an explanation would be of great help. Thank you for your time.

Answer 1

Spectral methods are global:

• Modifying a single value automatically changes the whole approximation

If you force the Dirichlet BC on the left side you unfortunately change the derivative (Neumann BC) on the right side. This has to be corrected by calculating the non-zero gradient on the right side and afterwards subtracting the weighted value from the right DOF.

The inverse weighting is motivated by the fact, that the gradient is calculated by a sum of collocation values times derivative matrix, globally. The last entry of $\mathbf{D}$ (last column, last line) is multiplied with the last value of $T$, which is responsible for the derivative at the last degree of freedom.

There should nothing change for the inhomogeneous case.

See Link