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.

  • 1
    $\begingroup$ The answer was given in another thread some days ago. scicomp.stackexchange.com/questions/43544/… $\endgroup$
    – ConvexHull
    Dec 21, 2023 at 6:11
  • 1
    $\begingroup$ Thank you for the reply @ConvexHull. It seems to me that you took the derivative of the solution vector at the boundary point, then you divided it with the corner value of the differentiation matrix and then subtracted it from the boundary point of the solution vector. Can you please explain how you derived that expression? I don't understand why that boundary point adjustment works. Why is the subtraction from the boundary point done like that and why do the dividing with the corner element of the differentiation matrix? What would change if the boundary condition wasn't homogeneous? $\endgroup$
    – Amilox Lex
    Dec 21, 2023 at 8:31

1 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

  • 1
    $\begingroup$ Thank you a lot for your explanation. I understand now. I substituted the corrected boundary value from your link into the numerical derivative at the last point and noticed that the result was zero. Suddenly, the correction step made sense :D $\endgroup$
    – Amilox Lex
    Dec 21, 2023 at 18:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.